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1 Introduction 

1.1 1925: Einstein's prediction for the ideal Bose gas 

Einstein considered N non-interacting bosonic and non-relativistic particles in a cubic 
box of volume L 3 with periodic boundary conditions. In the thermodynamic limit, defined 
as 

N 

N, L — > oo with — = p = constant, (1) 

Li 

a phase transition occurs at a temperature T c defined by: 

pA2 B (T c ) = C(3/2) = 2.612... (2) 

where we have defined the thermal de Broglie wavelength of the gas as function of the 
temperature T : 

and where ((a) = S&Li ^/k a is the Riemann Zeta function. 

The order parameter of this phase transition is the fraction N /N of particles in 
the ground state of the box, that is in the plane wave with momentum p = . For 
temperatures lower than T c this fraction N /N remains finite at the thermodynamic 
limit, whereas it tends to zero when T > T c : 

Nn 

T>T C ^0 (4) 

N _ /T\ 3 / 2 



For T <T C the system has formed a Bose-Einstein condensate in p = . The number 
iY of particles in the condensate is on the order of N , that is macroscopic. As we will 
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see, the macroscopic population of a single quantum state is the key feature of a Bose- 
Einstein condensate, and gives rise to interesting properties, e.g. coherence (as for the 
laser) . 

1.2 Experimental proof? 

The major problem encountered experimentally to verify Einstein's predictions is that at 
densities and temperatures required by Eq.(2) at thermodynamic equilibrium almost all 
materials are in the solid state. 

An exception is He 4 which is a fluid at T = . However He 4 is a strongly interacting 
system. In He 4 in sharp contrast with the prediction for the ideal gas Eq.(5), N Q /N < 
10% even at zero temperature [1]. 1 

The solution which victoriously led to Bose-Einstein condensation in atomic gases is 
to bring the system to extremely low densities (much lower than in a normal gas) and to 
cool it rapidly enough so that it has no time to recombine and solidify. The price to pay 
for an ultralow density is the necessity to cool at extremely low temperatures. Typically 
one has in the experiments with condensates: 



The critical temperatures range from 20 nK to the \i K range. 

Bose-Einstein condensation was achieved for the first time in atomic gases in 1995. 
The group of Eric Cornell and Carl Wieman at JILA was first, with 87 Rb atoms [2]. 
They were closely followed by the group of Wolfgang Ketterle at MIT with 23 Na atoms 
[3] and the group of Randy Hulet at Rice University with 7 Li atoms [4] . Nowadays there 
are many condensates mainly with rubidium or sodium atoms. No other alkali atoms 
than the ones of year 1995 has been condensed. Atomic hydrogen has been condensed in 
1998 at MIT in the group of Dan Kleppner [5]; the experiments on hydrogen were actually 
the first ones to start and played a fundamental pioneering role in developing many of 

1 Amusingly the ideal gas prediction Eq.(2) does not give a too wrong result for the transition temper- 
ature in helium. Note that the condensate fraction Nq/N should not be confused with the superfluid 
fraction: at T = the superfluid fraction is equal to unity. 



p < 10 15 atoms/cm 3 
T < lfiK. 



(6) 
(?) 
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the experimental techniques having led the alkalis atoms to success, such as magnetic 
trapping and evaporative cooling of atoms. 

In our lectures we do not consider the experimental techniques used to obtained and to 
study Bose-Einstein condensates as they are treated in the lectures of Wolgang Ketterle 
and of Eric Cornell at this school. 

1.3 Why interesting? 

1.3.1 Simple systems for the theory 

An important theoretical frame for Bose-Einstein condensation in interacting systems 
was developed in the 50's by Beliaev, Bogoliubov, Gross, Pitaevskii in the context of 
superfluid helium. This theory however is supposed to work better if applied to Bose 
condensed gases where the interactions are much weaker. 

The interactions in ultracold atomic gases can be described by a single parameter 
a , the so-called scattering length, as interactions take place between atoms with very 
low relative kinetic energy. The gaseous condensates are dilute systems as the mean 
interparticle separation is much larger than the scattering length a : 

p|a| 3 <l. (8) 

This provides a small parameter to the theory and, as we shall see, simple mean field 
approaches can be used with success to describe most of the properties of the atomic 
condensates. 

1.3.2 New features 

Atomic gases offer some new interesting features with respect to superfluid helium 4: 

• Spatial inhomogeneity: This feature can be used as a tool to detect the presence of 
a Bose-Einstein condensate inside the trap: in an inhomogeneous gas Bose-Einstein 
condensation occurs not only in momentum space but also in position space! 

• Finite size effects: The number of atoms in condensates of alkali gases is usually 
No < 10 7 . The hydrogen condensate obtained at MIT by Kleppner is larger 
N Q — 10 9 . Interesting finite size effects, that is effects which disappear at the 
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thermodynamic limit, such as Bose-Einstein condensates with effective attractive 
interactions ( a < ), can be studied in relatively small condensates. 

It is also interesting to consider small condensates where some interesting quan- 
tum aspects concerning coherence properties of the condensates, such as collapses 
and revivals of the relative phase between two condensates [6], could perhaps be 
measured [7]. 

• Tunability: Condensates in atomic gases can be manipulated and studied using the 
powerful techniques of atomic physics (see the lectures of Wolfgang Ketterle and 
Eric Cornell). Almost all the parameters can be controlled at will, including the 
interaction strength a between the particles. The atoms can be imaged not only 
in position space, but also in momentum space, allowing one to see the momentum 
distribution of atoms in the condensate! One can also tailor the shape and intensity 
of the trapping potential containing the condensate. 



2 The ideal Bose gas in a trap 

Let us consider a gas of non-interacting bosonic particles trapped in a potential U(f) 
at thermal equilibrium. As the particles do not interact thermal equilibrium has to be 
provided by coupling to an external reservoir. In the grand-canonical ensemble the state 
of the gas is described by the equilibrium N -body density matrix 

p=|exp[-/3(ff-MiV)] (9) 

where 5 is a normalization factor, H is the Hamiltonian containing the kinetic energy 
and trapping potential energy of all the particles, N is the operator giving the total 
number of particles, /3 = l/k B T where T is the temperature, and \i is the chemical 
potential. One more conveniently introduces the fugacity: 

z = exp[£/i]. (10) 

2.1 Bose-Einstein condensation in a harmonic trap 

Let us consider the case of a harmonic trapping potential U(f) : 

V{T) = \m{ufr?+<4tf + (11) 
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We wish to determine the properties of the trapped gas at thermal equilibrium; the 
calculations can be done in the basis of harmonic levels or in position space. 

2.1.1 In the basis of harmonic levels 

Let us consider the single particle eigenstates of the harmonic potential with eigenvalues 
€f labeled by the vector: 

I = {lx, l y , h) L = 0, 1, 2, 3... (a = x, y, z). (12) 

One has: 

€f= lxhW X + lyflUJy + l z huj z (13) 

where the zero-point energy (K/2)(oj x + u y + u z ) has been absorbed for convenience in 
the definition of the chemical potential. Let us consider the case of an isotropic potential 
for which all the uj a 's are equal to u , so that €f = Itwj with I = l x + l y + l z . 

The mean occupation number of each single particle eigenstate in the trap is given by 
the Bose distribution: 



1 



1 



i -i 



- exp(plhu) - 1 . (14) 



exp[£(e r -fj)]~ 1 Lz 

Since nj has to remain positive (for Z =0,1,2 ...), the range of variation of the fugacity 
z is given by 

0<z<l. (15) 

The average total number of particles N is obtained by summing over all the occupation 
numbers: N = Z)f^f > a relation that can be used in principle to eliminate z in terms 
of TV . It is useful to keep in mind that for a fixed temperature T , N is an increasing 
function of z . 

In the limit z — y one recovers Boltzmann statistics: rif cx exp(-/3e^) . We are 
interested here in the opposite, quantum degenerate limit where the occupation number 
of the ground state I = of the trap, given by 

No = n$ = (16) 
x z 

diverges when z — > 1 , which indicates the presence of a Bose-Einstein condensate in the 
ground state of the trap. 
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We wish to watch the formation of the condensate when z is getting closer to one, 
that is when one gradually increases the total number of particles N . The essence of 
Bose-Einstein condensation is actually the phenomenon of saturation of the population 
of the excited levels in the trap, a direct consequence of the Bose distribution function. 
Consider indeed the sum of the occupation numbers of the single particle excited states 
in the trap: 

JV = £«f- (17) 

The key point is that for a given temperature T , N' is bounded from above: 

N ' = H [I «p(#M - 1] 1 < E [exp(/3/M - i]- 1 = jv; ax . (is) 

Note that we can safely set z = 1 since the above sum excludes the term 1 = 0. 

If the temperature T is fixed and we start adding particles to the system, particles 
will be forced to pile up in the ground state of the trap when N > N max , where they will 
form a condensate. Let us now estimate the "critical" value of particle number AT max . 

We will restrict to the interesting regime k B T > Huj : in this regime Bose statistics 
allows one to accumulate most of the particles in a single quantum state of the trap 
while having the system in contact with a thermostat at a temperature much higher 
than the quantum of oscillation huj , a very counter-intuitive result for someone used 
to Boltzmann statistics! On the contrary the regime k B T <C huj would lead to a large 
occupation number of the ground state of the trap even for Boltzmann statistics. 

A first way to calculate N max is to realize that the generic term of the sum varies 
slowly with I as k B T ^> huj so that one can replace the discrete sum X)f^o by an 
integral J la>Q d 3 l . As we are in the case of a three-dimensional harmonic trap there is no 
divergence of the integral around 1 = 0. 

We will rather use a second method, which allows one to calculate also the first cor- 
rection to the leading term in ksT/huj . We use the series expansion 

1 p -x oo 

- x = Ee~ kX (19) 



e x - 1 1 - e~ x k=1 



which leads to the following expression for N max , if one exchanges the summations over 
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/ and k : 



N' 

max 



k=l i^q ol k=l 



1 — exp[— f3%ujk] t 



- 1 



(20) 



We now expand the expression inside the brackets for small x : 

3 

I -1 



1 



1 — exp[— x 
and we sum term by term to obtain 

NL„ = 



1 



x 6 



2x 2 



+ ... 



(21) 



(22) 



Note that the exchange of summation over k and summation over the order of expansion 
in Eq.(21) is no longer allowed for the next term 1/x , which would lead to a logarithmic 
divergence (that one can cut "by hand" at k ~ k B T/huj ). 

One then finds to leading order for the fraction of population in the single particle 
ground state: 



N~ N - 1 aS) [huj N ~ 1 



where the critical temperature T c ° is defined by: 

C(3) 



= N 



(23) 



(24) 



and C(3) = 1.202... . Note that the universal law (23) differs from the one obtained in 
the homogeneous case (5) usually considered in the literature. 

The present calculation is easily extended to the case of an anisotropic harmonic trap. 
To leading order one finds 



N r ~ 



k B T" 



C(3) 



(25) 



where uj = {uj x uj y uj z ) l / z is the geometric mean of the trap frequencies. One can also 
calculate A^ ax in two-dimensional and one-dimensional models. One also finds in these 
cases a finite value for N max : the saturation of population in the single particle excited 
states applies as well and one can form a condensate, a situation very different from the 
thermodynamical limit in the homogeneous ID and 2D cases. 
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2.1.2 Comparison with the exact calculation 

One can see in figure 1 that the first two terms in the expansion Eq.(22), combined with 
the approximation N Q /N ~ 1 - N^/N , give a very good approximation to the exact 
condensate fraction for N = 1000 particles only. 




T/T c ° 

Figure 1 : Condensate fraction versus temperature for an ideal Bose gas in a spherically symmet- 
ric trap with N = 1000 particles. The circles correspond to the exact quantum calculation. The 
solid line corresponds to the prediction Nq/N ~ with N^ ax given by the two terms 

in the expansion Eq.(22). The dashed line corresponds to the prediction Nq/N ~ 1 — N^^/N 
with AT^ax given by the leading term in Eq.(22). This figure was taken from [8]. 

2.1.3 In position space 

A very important object in the description of the state of the gas is the so-called one-body 
density matrix. We can define it as follows. 
Consider a one-body observable 



* = £*(») 

i=l 



(26) 
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where X(i) is the observable for particle number i and where N is the operator giv- 
ing the total number of particles. The one-body density matrix pi is defined by the 
requirement that for any X : 

M^TtbMTU)]. (27) 

For the particular case of X equal to the identity it follows X = N and (X) = Tr[/?i] = 
(N) so that our one-body density matrix is normalized to the mean number of particles 
in the system. 

An equivalent definition of p\ in the second quantized formalism is simply 

(f VilO = {^(f)W')) (28) 

where ^(r ) is the atomic field operator, annihilating an atom in r . 

At thermal equilibrium in the grand-canonical ensemble, the one-body density matrix 
of the ideal Bose gas is given by 

Pi = (29) 

z- 1 expiPhx) - 1 

where the single-particle Hamiltonian in the case of a spherically symmetric harmonic 
trap is 

hi = ^ h -muj 2 f 2 — -huj. (30) 

Here again we have subtracted the zero-point energy for convenience. The Bose formula 
Eq.(14) corresponds to the diagonal element of pi in the eigenbasis of the harmonic 
oscillator (the off-diagonal elements of course vanish). In position space the diagonal 
term 

{r\pi\r) = p(f) (31) 

gives the mean spatial density of the gas. 

In order to calculate the density we use the series expansion Eq.(19) to rewrite pi as 
follows: 

oo 

ft = £zV». (32) 

k=l 

This writing takes advantage of the fact that the matrix elements 

(r\e-f> khl \r') (33) 
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are known for an harmonic oscillator potential [9]. One then obtains explicitly: 



pin = ( 



^ 3/2 fy (1 - exp(-2^M)" 3/2 exp 



k=l 



mujr 



n 



tanh 



(34) 



One can identify the contribution of the condensate to this sum when z — > 1~ . When 
the summation index h is large, what determines the convergence of the series is indeed 
the factor z k . Replacing the other factors in the summand by their asymptotic value for 
k — ¥ +oo we identify the diverging part when z = 1 : 



) z ex P 



nth 



k=l 



mujr 



1 -z 



|0oao(O| 2 = Wo A o(Of 



(35) 



where 0o,o,o(O is the ground state wave function of the harmonic oscillator. 

Numerically we have calculated the total density p(f) for a fixed temperature k B T = 
20huj and for increasing number of particles (see figure 2). Here the maximal number of 
particles one can put in the excited states of the trap is ~ ((3)(kBT/hio) 3 ~ 10 4 . 

When N <§C Nmax the effect of an increase of AT is mainly to multiply the density 
by some global factor (the curves in logarithmic scale in figure 2 are parallel one to 
the other). When TV is becoming larger than A^ ax a peak in density grows around 
r = , indicating the formation of the condensate, whereas the far wings of the density 
distribution saturate, which reflects the saturation of the population of the excited levels 
of the trap. 



2.1.4 Relation to Einstein's condition p\ z dB = ((3/2) 

In the limit k B T ^> Tiuj we can actually calculate the value p f max (r ) to which the density 
p'(f) of particles in the excited states of the trap saturates when z — » 1 . We simply use 
the expansion Eq.(34), subtracting from the total density p(f ) the contribution of the 
condensate iV o |0o,o,o(^)| 2 ■ The resulting series is converging even for z = 1 so that we 
can take safely the semiclassical limit k B T ^> Tiuj term by term in the sum: 



p'(f ) ~ -3 



E 



a5b *f i * 3/2 



exp 



zexp ( —d-muPr 2 



where 



oo x k 



(36) 



(37) 



k=l 
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10 4 




r[aj 



Figure 2: Spatial density for an ideal Bose gas at thermal equilibrium in a harmonic trap of 
frequency uj . The temperature is fixed to ksT = 20huj and the number of particles ranges 
from N = 500 to N = 32000 between the lowest curve and the upper curve, with a geometrical 
reason equal to 2. The unit of length for the figure is ao = {h/2mu>) 1 / 2 , that is the spatial 
radius of the ground state of the trap. 

We term this approximation semiclassical as (i) one can imagine that the classical limit 
ft — > is taken in each term k of the sum, giving the usual Gaussian distribution for the 
density of a classical harmonic oscillator at temperature ksT/k , but (ii) the distribution 
still reflects the quantum Bose statistics. 

If now we set z = 1 in (36) to express the fact that a condensate is formed we obtain 

p' m ^(f= 0) ~ ^ 3/2 (l) = (38) 

A dB A dB 

We therefore recover Einstein's condition provided one replaces the density p of the 
homogeneous case by the density at the center of the trap. 

2.2 Bose-Einstein condensation in a more general trap 

We now extend the idea of the previous semiclassical limit to more general non-harmonic 
potentials. This allows to find the condition for Bose-Einstein condensation in presence 
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of a non-harmonic potential. This will prove useful in presence of interactions between 
the particles where the non-harmonicity is provided by the mean field potential. 



2.2.1 The Wigner distribution 

The idea is to find a representation of the one-body density matrix having a simple (non 
pathological) behavior when H — »• . Let us take as an example a single harmonic 
oscillator. The density matrix is then of the form: 



a = — e~ 0H ^o 
Z 



(39) 



where H ho is the harmonic oscillator Hamiltonian. As shown in [9] all the matrix elements 
of a can be calculated exactly: 



(flair') = 



1 



(27r) 3 / 2 (Ar) 3 



exp 



[(r + f')/2F 



2(Ar) 5 





" (f_f')2- 


exp 


[ 2 £ 2 J 



The relevant length scales are the spatial width of the cloud Ar : 



(Ar) 2 = — — cotanh 
v 1 2muj 



2k B T 



and the coherence length f : 



(40) 



(41) 



(42) 



If we now take the classical limit h -> (in more physical terms the limit hu < k B T ') 
then: 

(Ar) 2 M (43) 

h 2 \ 2 

e ~ _J^ = ^s. (44) 

^ mk B T 2tt v ; 

In the limit ft — >• the ft dependence of £ causes (f \a\f ') — > for fixed values of r, f' 
unless f=r f : the limit is singular. 

To avoid this problem one can use the Wigner representation of the density matrix, 
introduced also in the lectures of Zurek and Paz: 

d 3 u 



W[a](r,p) = J 



a 



2 



(45) 
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The Wigner distribution is the quantum analog of the classical phase space distribution. 
In particular one can check that the Wigner distribution is normalized to unity and that 



j d 3 fW(f,p) = (p\°\p) 
J d 3 pW(r,p) = (f \<r\r). 



(46) 
(47) 



An important caveat is that W is not necessarily positive. 

For the harmonic oscillator at thermal equilibrium the integral over u in Eq.(45) is 
Gaussian and can be performed exactly: 



W(f,p) = 



1 



exp(- 



(2^ArApf ^ 2(Ar) 2 
where Ap = h/£ . If we take now the limit % — » : 



) exp(- 



r 



2(A P y 







(Arf 



moo* 



(48) 

(49) 
(50) 



(Ap) 2 -> mk B T 
so that W(f,p) tends to the classical phase space density. 

2.2.2 Critical temperature in the semiclassical limit 

Let us turn back to our problem of trapped atoms in a non-harmonic trap where the 
single particle Hamiltonian is given by 

k = ^ + U(f) (51) 
and the one-body density matrix is given by Eq.(32). For ^40 we have: 

- k p(£ + U(?))}. (52) 



W[e- pkhl ](f,p) ~ ^exp 



As we did before we put apart the contribution of the condensate. One then gets for 
the one-body density matrix of the non-condensed fraction of the gas in the semiclassical 
limit: 



1 f 1 



/i 3 



exp 



-i 



- 1 



(53) 
(54) 
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We are now interested in the spatial density of the non-condensed particles in the 
semiclassical limit. By integrating Eq.(53) over p we obtain: 

P'sdr) = v^3/ 2 (ze-W>) (55) 

A dB 

where g a is defined in Eq.(37). The condition for Bose-Einstein condensation is z ^ 
e^^min where U m i n = min^ U(f ) is the minimal value of the trapping potential, achieved 
in the point r min . For z = e^min the semiclassical approximation for the non-condensed 
density gives in this point: 

/4(r min ) = -73-03/2(1) (56) 

A dB 

or 

pX 3 dB = 2.612... (57) 

Again Einstein's formula is recovered with p being the maximal density of the non- 
condensed cloud, that is the non-condensed density at the center of the trap. 

The semiclassical calculation that we have just presented was initially put forward 
in [10]. We do not discuss in details the validity of this semi-classical approximation. 
Intuitively a necessary condition is ksT ^> AE where AE is the maximal level spacing 
of the single particle Hamiltonian among the states thermally populated. Some situations, 
where the trapping potential is not just a single well, may actually require more care. The 
case of Bose-Einstein condensation in a periodic potential is an interesting example that 
we leave as an exercise to the reader. 

2.3 Is the ideal Bose gas model sufficient: experimental verdict 
2.3.1 Condensed fraction as a function of temperature 

The groups at MIT and JILA have measured the condensate fraction N Q /N as function 
of temperature for a typical number of particles N — 10 5 or larger. We reproduce 
here the results of JILA [11] (see figure 3). This figure shows that the leading order 
prediction of the ideal Bose gas Eq.(23) is quite good, even if there is a clear indication 
from the experimental data that the actual transition temperature is lower than T c ° . This 
deviation may be due to finite size effects and interaction effects but the large experimental 
error has not allowed yet a fully quantitative comparison to theory. 
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Figure 3: Condensate fraction N Q /N as function of T/T c ° where T c ° is the leading order 
ideal Bose gas prediction Eq.(24). Circles are the experimental results of [11] while the dashed 
line is Eq.(23). 

2.3.2 Energy of the gas as function of temperature and number of particles 

In the experiments one produces first a Bose condensed gas at thermal equilibrium. Then 
one switches off suddenly the trapping potential. The cloud then expands ballistically, 
and after a time long enough that the expansion velocity has reached a steady state value 
one measures the kinetic energy of the expanding cloud. 

Suppose that the trap is switched off at t = . For t = 0~ the total energy of the 
gas can be written as 

^tot(0~) — ^kin + ^trap + ^intj (58) 

that is as the sum of kinetic energy, trapping potential energy and interaction energy. At 
time t = + there is no trapping potential anymore so that the total energy of the gas 
reduces to 

^tot(0 + )=^ kin + ^inf (59) 

In the limit t —¥ +00 the gas expands, the density and therefore the interaction energy 
drop, and all the energy £^ t ot(0 + ) is converted into kinetic energy, which is measured. 
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In figure 4 we show the results of JILA for E tot (0 + ) for temperatures around T c ° 
[11] together with the ideal Bose gas prediction. The main feature of the ideal Bose gas 
prediction is a change in the slope of the energy as function of temperature when T 
crosses T c . One observes indeed a change of slope in the experimental results (see the 
magnified inset)! 

For T > T c the ideal Bose gas model is in good agreement with the experiment. For 
T <T C we observe however that the experiment significantly deviates from the ideal Bose 
gas. 




Figure 4: Expansion energy of the gas E tot (0 + ) per particle and in units of ksT® as function 
of the temperature in units of T c ° . The disks correspond to the experimental results of [11]. 
The straight solid line is the prediction of Boltzmann statistics. The dashed curve exhibiting a 
change of slope is the ideal Bose gas prediction. The curved solid line is a piecewise polynomial 
fit to the data. The inset is a magnification showing the change of slope of the energy as function 
of T close to T = T c ° . The figure is taken from [11]. 

What happens at even lower values of T/T c ° ? We show in figure 5 the expansion 
energy of the condensate per particle in the regime of an almost pure condensate [12]. 
This energy then depends almost only on the number of condensate particles N Q , in a 
non-linear fashion. This is in complete violation with the ideal Bose gas model, which 
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predicts an energy per particle in the condensate independent of N Q . More precisely 
the ideal Bose gas prediction would be h(u x + uj y + oj z )/4: where the uj a 's are the trap 
frequencies. In units of ks this would be in the 10 nK range, an order of magnitude 
smaller than the measured values. 
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Figure 5: Expansion energy of the condensate per particle in the condensate, divided by ks , 
as a function of the number of particles in the condensate. The experiment is performed at 
temperatures T « T c . The triangles correspond to cases where the non-condensed cloud was 
not visible experimentally. The disks correspond to case where the non-condensed cloud could 
be seen. The figure is taken from [12]. The solid line is a fit of the interacting Bose gas prediction 
of §5. 

2.3.3 Density profile of the condensate 

The group of Lene Hau at Rowland Institute has measured the density profile of the 
condensate in a cigar-shaped trap, along the weakly confining axis z of the trap. As 
imaging with a light beam is used the actual density obtained in the experiment is the 
density integrated along the direction y of propagation of the laser beam, plotted in 
figure 6 for x = as function of z [13]. The measured profile is very different from and 
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much broader than the Gaussian density profile of the ground state wavefunction of the 
harmonic oscillator. 
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Figure 6: Column density profile (see text) of a condensate along the weak axis z of a cigar- 
shaped trap. The experimental results of [13] (dots) are very different from the ideal Bose gas 
prediction (dashed line). The solid line corresponds to the theoretical prediction of §5. 

2.3.4 Response frequencies of the condensate 

By modulating the harmonic frequencies of the trapping potential one can excite breathing 
modes of the condensate. For example the group at MIT modulated the trap frequency 
along the slow axis z of a cigar-shaped trap and observed at T <C T c subsequent 
breathing of the condensate at a frequency 1.569 (4)0^ . This frequency is not an integer 
multiple of u z and can therefore not be obtained in the ideal Bose gas model. 

In conclusion the ideal Bose gas model may be acceptable as long as no significant 
condensate has been formed. If a condensate is formed interaction effects become impor- 
tant, and dominant at T «C T c . This serves as a motivation to the next sections of this 
lecture, which will deal with the interacting Bose gas problem. 
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3 A model for the atomic interactions 

The previous section 2 has shown that the ideal Bose gas model is insufficient to explain 
the experimental results when a condensate is formed. In this section we choose the 
model potential to be used in this lecture to take into account the atomic interactions. 
The reader interested in a more careful discussion of real interaction potentials is referred 
to [14]. 

3.1 Reminder of scattering theory 

We consider two particles of mass m interacting in free space via the potential V(f[ — f 2 ) 
depending on the positions f[ , rj only through the relative vector f[ — ■ The center of 
mass of the two particles is then decoupled from their relative motion, and the evolution 
of the relative motion is governed by the Hamiltonian: 

H ^ = ^ + V(f) (60) 

where r = f[ — r~2 is the vector of coordinates of the relative motion, p = (p[ — p^)/2 is 
the relative momentum and \i — m/2 is the reduced mass. We assume in what follows 
that the potential V(f) is vanishing in the limit r — »• oo . 

3.1.1 General results of scattering theory 

The scattering states ip{r) of the relative motion of the two particles are the eigenstates 
of H re \ with positive energy E . Writing E — h 2 k 2 /2}i and multiplying the eigenvalue 
equation by 2fi/h 2 we obtain 

(A + k^{r) = 2 ^V{r)i>{r). (61) 

One has also to specify boundary conditions on tp to get the full description of a scattering 
state. This is achieved by means of an integral formulation of the eigenvalue equation. 

• Integral equation 

To obtain the integral formulation of the scattering problem we write the right hand side 
of the eigenvalue equation Eq.(61) as a continuous sum of Dirac distributions: 

(A + k 2 )i/j(r ) = J dV'pV(?0^W- ?')- ( 62 ) 
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We then find a solution of this equation with a single Dirac distribution on the right hand 
side: 

{^+k 2 )i) G {?) = (63) 
having the form of an outgoing spherical wave for r — > oo : 

i p ik\f— f'\ 

Mr) = -2-T"7P ( 64 ) 

This is actually a Green's function of the operator A + k 2 . The scattering state of the 
full problem can then be written as 

9// r p ik\f-f'\ 

m = Mr) - ^2 / d 3 r'y-^{? ')W •) . (65) 

The first term tp is the incoming free wave of the collision, solving (A + k 2 )tp Q = ; we 

— # 

simply assume here that the incoming wave is a plane wave of wavevector k : 

ipo(r) = exp[ik-r]. (66) 
The remaining part of ip is then simply the scattered wave. 

• Born expansion 

When the interaction potential is weak one sometimes expands the scattering state ip 
in powers of V . In the integral formulation Eq.(65) of the eigenvalue equation this 
corresponds to successive iterations of the integral, the approximation for ip at order 
n + 1 in V being obtained by replacing ip by its approximation at order n in the 
right-hand side of the integral equation. E.g. to zeroth order in V , ip — ip Q , and to first 
order in V we get the so-called Born approximation: 

9// r p ik\f-f'\ 

V'Bom {r) = Mr) -Tit \ d 3 r' ^^V(f')Mr'). (67) 

AtyH J \r — r'\ 

3.1.2 Low energy limit for scattering by a finite range potential 

Some results can be obtained in a simple way when the potential V has a finite range 
b , that is when it vanishes when r > b . 

• asymptotic behavior for large r 
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As the integration over the variable r ' is limited to a range of radius b one can expand 
the distance from r to f ' in powers of r when r > b : 

\f-f'\ = r -7 f .fi + (68) 

where ft = r/r is the direction of scattering. The neglected term, scaling as b 2 /r , has a 
negligible contribution to the phase exp[ik\f - f'\] when r > kb 2 . One then enters the 
asymptotic regime for ip : 

m = Mr) + + O (^) (69) 

where the factor , the so-called scattering amplitude, does not depend on the distance 
r : 

h[n) = / d¥ ' e ~ im ' f ' v (r'mn (70) 

If the mean distance between the particles in the gas, on the order of p -1 / 3 , where 
p is the density, lies in the asymptotic regime for ip (that is p~ 1 ^ > b, kb 2 ) the effect 
of binary interactions on the macroscopic properties of the gas will be sensitive to the 
scattering amplitude f% , and no longer to the details of the scattering potential. This 
is the key property that we shall use later in this low density regime to replace the 
exact interaction potential by a model potential having approximately the same scattering 
amplitude. 

• limit of low energy collisions 

Another simplification comes from the fact that collisions take place at low energy in 
the Bose condensed gases: as f?k 2 is on the order of k B T in the thermal gas, k 
becomes small at low temperature. 

If kb<^l the phase factor exp[—ikn-f '] becomes close to one in the integral Eq. (70) 
giving the scattering amplitude. The scattering amplitude then no longer depends on 
the scattering direction ft , the asymptotic part of the scattered wave becomes spherically 
symmetric (even if the scattering potential is not!): one then says that scattering takes 
place in the s -wave only. 

Going to the mathematical limit k — > we get for the scattering amplitude: 



(71) 
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The quantity a is the so-called scattering amplitude; it will be the only parameter of our 
theory describing the interactions between the particles, and our model potential will be 
adjusted to have the same scattering length as the exact potential. When k is going to 
zero, the scattering state converges to the zero energy scattering state, behaving for large 
r as 



A numerical calculation of this zero energy scattering state is an efficient way of calculating 
a for a given potential V . Note that there is of course no connection between a and 
b , except for particular potentials like the hard sphere potential. 

3.1.3 Power law potentials 

In real life the interaction potential between atoms is not of finite range, as it contains 
the Van der Waals tail scaling as 1/r 6 for large r 2 . It is fortunately possible to show 
for the class of power-law potentials, scaling as l/r n , that several of our conclusions, 
obtained in the finite range case, hold provided that n > 3 . E.g. in the limit of small 
k 's only the s -wave scattering survives, and fj* has a well defined limit for k — > , 
allowing one to define the scattering length. 

3.2 The model potential used in this lecture 
3.2.1 Why not keep the exact interaction potential ? 

For alkali atoms the exact interaction potential has a repulsive hard core, is very deep (as 
deep as 10 3 Kelvins times ks for 133 Cs), has a minimum at a distance r i2 on the order 
of 6 A(for cesium), and contains many bound states corresponding to molecular states of 
two alkali atoms (see figure 7). 

There are several disadvantages to use the exact interaction potential in a theoretical 
treatment of Bose-Einstein condensation: 

1. V is difficult to calculate precisely, and a small error on V may result in a large 
error on the scattering length a . In practice a is measured experimentally, and 
this is the most relevant information on V in the low density, low temperature 
limit. 

2 or even as 1/r 7 if r is larger than the optical wavelength. 




(72) 



Model for interactions 



27 




-750 1 ■ 1 ' 1 ■ 1 

0.5 1 1.5 

r 12 [nm] 

Figure 7: Typical shape of the interaction potential between two atoms, as function of the 
interatomic distance Tyi . The numbers are indicative and correspond to cesium. 

2. the presence of bound states of V with a binding energy much smaller than the 
temperature of the gas (there are 9 orders of magnitude between the potential depth 
10 3 K and the gas temperature ~ 1/i K) clearly indicates that the Bose condensed 
gases are in a metastable state; at the experimental temperatures and densities 
the complete thermal equilibrium of the system would be a solid. Direct thermal 
equilibrium theory, such as the thermal N -body density matrix exp[— /3H] , cannot 
therefore be used with V . This is why even in the exact Quantum Monte Carlo 
calculations performed for alkali gases [15] V is replaced by a hard sphere potential. 
Such a complication was absent for liquid helium, where the well-known exact V 
can be used [16]. 

3. V can not be treated in the Born approximation, because it is very strongly re- 
pulsive at short distances and has many bound states: even if the scattering length 
was zero, one would have to resum the whole Born series to obtain the correct result 
[We recall that for a potential as gentle as a square well of radius b , the Born ap- 
proximation applies when the zero-point energy for confinement within a domain of 
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radius b , h 2 /2fib 2 , is much larger than the potential depth, which implies that no 
bound state is present in the well.] As a consequence naive mean field approxima- 
tions, which neglect the correlations between particles due to interactions, implicitly 
relying on the Born approximation, cannot be used with the exact V . 

The key idea is therefore to replace the exact interaction potential by a model potential 
(i) having the same scattering properties at low energy, that is the same scattering length, 
and (ii) which should be treatable in the Born approximation, so that naive mean field 
approaches apply. 

The model potential satisfying these requirements with the minimal number of pa- 
rameters (one!) is the zero-range pseudo-potential initially introduced by Enrico Fermi 
[17, 18] and having the following action on any two-body wavefunction: 



(fl,fi\V\ipi } 2) = gS(f[ - rj) 



_dru 

The factor g is the so-called coupling constant 



(73) 

J ri2=0 



g = a 74 

m 

where a is the scattering length of the exact potential. The pseudo-potential involves a 
Dirac distribution and a regularizing operator. 

• Effect of regularization 

When the wavefunction ^ 1>2 is regular close to f[ = f % , one can check that the regular- 
izing operator has no effect, so that the pseudo-potential can be viewed as a mere contact 
potential g5(f[ — rj) . 

When the wavefunction ^ 1>2 has a l/ri2 divergence: 

Wi,2 (ri, r 2 ) = h regular (75) 

r 12 

where A is the function of the center of mass coordinates only the regularizing operator 
removes the diverging part: 

°(r a «*m=0. (76) 
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In this way we have extended the Hilbert space of the state vectors of the particles 
with wave functions diverging as l/ri 2 ; note that these wavefunctions remain square 
integrable, as the element of volume scales as rf 2 in 3D. As we shall see this l/ru 
divergence is a consequence of the zero-range of the pseudo-potential. 



3.2.2 Scattering states of the pseudo-potential 



Turning back to the relative motion of two particles we now derive the scattering states of 
the pseudo-potential from the integral equation Eq.(65). As the pseudo-potential involves 
a Dirac S(f f ) the integral over f' can be performed explicitly: 



d 



dr 



(rV(f ')) 



Jr'=0 



As the factor 



C = 



d_ 

dr 



7 (»W)) 



(77) 



(78) 



J r'=0 



does not depend on r we find that ip has the standard asymptotic behavior of a scat- 
tering state in r but everywhere in space, not only for large r . This is due to the 
zero-range of the pseudo-potential. To calculate C , we multiply Eq.(77) by r , we take 
the derivative with respect to r and set r to zero. On the left hand side we recover the 
constant C by definition. We finally obtain: 



C=l- aCik 



(79) 



so that C = 1/(1 + ika) and the scattering states of the pseudo-potential are exactly 
given by 



fcV J l + ika r 

The corresponding scattering amplitude, 

fk = 



1 + ika 



(81) 



does not depend on the direction of scattering, so that the pseudo-potential scatters only 
in the s -wave, whatever the modulus k is. The scattering length of the pseudo-potential, 
—fk=o = a , coincides with the one of the exact potential. 
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Finally we note that the total cross-section for scattering of identical bosons by the 
pseudo-potential is given by a Lorentzian in k , 

and that the pseudo-potential obeys the optical theorem. 
3.2.3 Bound states of the pseudo-potential 

As a mathematical curiosity we now point out that not only the scattering states but also 
the bound states of the pseudo-potential can be calculated. A first way of obtaining the 
bound states is a direct solution of Schrodinger's equation. A more amusing way is to use 
the following closure relation: 

(f^l^XI^Nl-Pbound (83) 

where \tp^) is the scattering state given in Eq.(80) and Pbound is the projector on the 
bound states of the pseudo-potential. 

In calculating the matrix elements of this closure relation between perfectly localized 
state vectors \f ) and |f ') and using spherical coordinates for the integration over k 
one ultimately faces the following type of integrals: 

/+oo F ik(r+r') 
dk — . (84) 
-oo 1 + ika 

We calculate / using the residues formula, by extending the integration variable k to 
the complex plane and closing the contour of integration by a circle of infinite radius, 
which has to be in the upper half of the complex plane as r + r' > . As the integrand 
in / has a pole in k = i/a , we find that / vanishes for a < , as the pole is then in 
the lower half of the complex plane. For a > the pole gives a non-zero contribution to 
the integral: 

J = 2 JL e -(rW)la^ (g5) 
d 

Finally we find that Pbound = for a < , corresponding to the absence of bound 
states, and -Pbound — Inbound) (Abound | & T a > 0, corresponding to the existence of a 
single bound state: 

1 e _r / a 

Abound (r ) = —7= . (86) 



V2 



ma 
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From Schrodinger's equation, we find for the energy of the bound state: 

ft 2 

Abound — o- (87) 

ma 2 

The existence of a bound state for a > and its absence for a < is a paradoxical 
situation. As we shall see in the mean field approximation, the case a > corresponds to 
effective repulsive interactions between the atoms, whereas the case a < corresponds 
to effective attractive interactions. In the purely ID case, the situation is more intuitive, 
the potential gioS(x) having a bound state only in the effective attractive case gw < . 
This paradox in 3D comes from the non-intuitive effect of the regularizing operator (an 
operation not required in ID), which makes the pseudo-potential different from a delta 
potential; actually one can shown in 3D that a delta potential viewed as a limit of square 
well potentials with decreasing width b and constant area does not scattered in the limit 
&->0. 



3.3 Perturbative vs non-perturbative regimes for the pseudo- 
potential 

3.3.1 Regime of the Born approximation 

As we will use mean field approximations requiring that the scattering potential is treat- 
able in the Born approximation, we identify the regime of validity of the Born approxi- 
mation for the pseudo-potential. 

As we have seen in the previous subsection the integral equation for the scattering 
states of the pseudo-potential can be reduced to the equation for C : 

C=l- ikaC, (88) 

the scattering state being given by 

tj) % {f) = e^' ? -aC^. (89) 

The Born expansion will then reduces to iterations of Eq.(88). To zeroth order in the 
interaction potential, we obtain C = so that reduces to the incoming wave. To 
first order, we get the Born approximation 

C\ = l- ikaC Q = 1. (90) 
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To second order and third order we obtain 

C 2 = 1 - ikaC x = 1 - ika (91) 
63 = 1- ikaC 2 = 1 - ika + (i/ca) 2 (92) 

so that the Born expansion is a geometrical series expansion of the exact result C — 
1/(1 + ika) in powers of ika . 

The validity condition of the Born approximation is that the first order result is a 
small correction to the zeroth order result. For the scattering amplitude this requires 

k\a\ < 1. (93) 

For the scattering state this requires 

r > a. (94) 

If one takes for r the typical distance p~ 1 ^ between the particles in the gas, where p 
is the density, this leads to 

p 1/3 \a\ < 1. (95) 

• Are the conditions for the Born approximation satisfied in the experiments ? 

To estimate the order of magnitude of k we average k 2 over a Maxwell-Boltzmann 
distribution of atoms with a temperature T — lp K, typically larger than the critical 
temperature for alkali gases; the average gives a root mean square for k equal to 

M= (ir) ■ (96) 

For 23 Na atoms used at MIT, with a scattering length of 50 a Bo hr , where the Bohr radius 
is ciBohr = 0.53 A, we obtain AA; a = 2 x 10~ 2 . For rubidium 87 Rb atoms used at JILA, 
with a scattering length of 110a Bo hr , we obtain Ak a = 0.1 . 

In the case of an almost pure condensate in a trap, the typical k is given by the 
inverse of the size R of the condensate, as the condensate wavefunction is not very far 
from a minimum uncertainty state. Generally this results in a much smaller AA; than 
Eq.(96), as R is much larger than the thermal de Broglie wavelength. One could however 
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imagine a condensate in a very strongly confining trap, such that R would become close 
to a ; in this case, not yet realized, the mean field theory has to be revisited. 

We turn to the second condition Eq.(95). The typical densities of condensates are on 
the order of 2 x 10 14 atoms per cm 3 . For the scattering length of sodium this leads 
to p l ^a ~ 0.015 <C 1 . For the scattering length of rubidium this leads to p l ^a ~ 
0.034 <C 1 . Both conditions for the Born approximation applied to the pseudo-potential 
are therefore satisfied. 

3.3.2 Relevance of the pseudo-potential beyond the Born approximation 

Let us try to determine necessary validity conditions for the substitution of the exact 
interaction potential by the pseudo-potential. 

First one should be in a regime dominated by s -wave scattering, as the pseudo- 
potential neglects scattering in the other wave. This condition is easily satisfied in the 
\i K temperature range for Rb, Na. 

Second the scattering amplitude of the exact potential in s -wave should be well 
approximated by the pseudo-potential. For isotropic potentials vanishing for large r as 
\jr n , with n > 5 , the s -wave scattering amplitude has the following low k expansion: 

f *~° = ~ a ^ + ik-\k 2 r e + ... (97) 

where r e is the so-called effective range of the potential. To this order in k the result of 
the pseudo-potential corresponds to the approximation r e = . When r e is on the order 
of a (which is the case for a hard sphere potential, but not necessarily true for a more 
general potential) the term in r e can be neglected if k 2 r e < that is (ka) 2 <C 1 ; 
there is therefore no meaning to use the pseudo-potential beyond the Born regime. 

Consider now the case r e <C \a\ . The term r e k 2 remains small as compared to 1/a 
for k\a\ < 1 . For k\a\ ^> 1 the term ik dominates over 1/a ; k 2 r e remains small as 
compared to ik as long as kr e <§C 1 . The use of the pseudo-potential may then extend 
beyond the Born approximation. 

An example of a situation with r e <§C \a\ is the so-called zero energy resonance, where 
a is diverging. When a bound state of the interaction potential is arbitrarily close to the 
dissociation limit, the scattering length diverges a —¥ +oo , the bound state has a large 
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tail in r scaling as e r l a /r and the bound state energy scales as —% 2 /ma 2 [19, 20]. 
These scaling laws hold for the pseudo-potential, as we have seen. 

4 Interacting Bose gas in the Hartree-Fock approxi- 
mation 

Now that we have identified a simple model interaction potential treatable in the Born 
approximation we use it in the simplest possible mean field approximation, the so-called 
Hartree-Fock approximation. This approximation was applied to trapped gases for the 
first time in 1981 (see [21])! 



The Hartree-Fock mean field approximation can be implemented in a variety of ways. We 
have chosen here the approach in terms of the BBGKY hierarchy, truncated to first order. 

4.1.1 Few body- density matrices 

We have already introduced in §2 the concept of the one-body density matrix. We revisit 
here this notion and extend it to two-body density matrices. 

• For a fixed total number of particles 

Let us first consider a system with a fixed total number of particles N and let <7i ;2 ...iv 
be the N -body density matrix. Starting from gi^...n we introduce simpler objects as 
the one-body and two-body density matrices pi and pn , by taking the trace over the 
states of all the particles but one or two: 



In practice the knowledge of pi and p i2 is sufficient to describe most of the experimental 
results. As you know, (r|pi|r) is the density of particles and (fi,f* 2 \pi\f{,f*2) is the pair 
distribution function. 



4.1 BBGKY hierarchy 



p[ = iVTr 2 ,3... JV (ai ; 2,...jv) 

p£> = N(N-l)Tr 3 ^ N (ai, 2 ,... N ). 



(98) 
(99) 
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• For a fluctuating total number of particles 

If N fluctuates according to the probability distribution Pn , we define few-body density 
matrices by the following averages over N : 

pi = E p »p { i N) (ioo) 

N 

Pn = Y, P NP { 12 ] (101) 

N 

Alternatively on can define directly the one-body and two-body density matrices in 
second quantization: 

(ri\pi\r 2 ) = (^(^(fl)) (102) 
(f[,r 2 \pi2\r^,fl) = (^(^^(rl)^^)^^!)). (103) 

Note that the few-body density matrices are normalized as 

Tr[^] = (N) (104) 
Tr[/3 12 ] = {N(N-l)) (105) 

so that one can obtain the variance of the fluctuations in the number of atoms from the 
one-body and two-body density matrices. 

4.1.2 Equations of the hierarchy 

The idea of our derivation of the mean field approximation is to get an approximate closed 
equation for pi by closing the hierarchy with some "cooking recipe" giving /? 12 in terms 
of p x . 

To derive the first equation of the hierarchy we start from the exact master equation: 

ih^0i,2..N = [H, a h2 .. N ] (106) 
where the Hamiltonian is the sum of one-body and two-body terms: 

H = J2h i + 1 - E % (107) 

i=l i^j-ij = l 
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The single particle Hamiltonian hi contains the kinetic and trapping potential energy of 
the atom i and Vij in the interaction potential between the atoms i and j . Now we 
take the trace of the master equation over the particles 2,3 ... N and multiply it by N , 
obtaining 



We have kept here only the terms involving the atom 1, as the other terms are commutators 
of vanishing trace. The sum over j amounts to N — 1 times the same contribution, e.g. 
the j = 2 contribution, as the atoms are indiscernible. We finally obtain the first equation 
of the hierarchy: 



The equation Eq.(109) is not closed for pi , as it involves p i2 . The next equation of 
the hierarchy, the equation for p 12 , involves p i23 , etc, up to the N -body density matrix, 
where the hierarchy terminates. The mean field approximation consists in replacing p i2 
by an ad hoc function of pi . 

4.2 Hartree-Fock approximation for T > T c 
4.2.1 Mean field potential for the non-condensed particles 

We use the following simple approximation to break the hierarchy: 



where P i2 is the permutation operator exchanging the states of the particles 1 and 2. 
The last identity in (110) is obtained by using the commutation of F i2 and pi® pi , and 
the fact that Ff 2 = 1 . 

The factorized prescription p i2 — pi® pi is the Hartree approximation. It assumes 
weak correlations between the particles. Indeed at short distances r i2 , the real p i2 
is expected to be a statistical mixture of scattering states of the interaction potential. 
Neglecting the correlations in p 12 between particles 1 and 2 amounts to considering 
only separable, plane wave scattering states, which corresponds to the zeroth order in 
the Born expansion of the scattering theory. Actually p i2 appears in Eq.(109) inside a 
commutator with V 12 , so that taking the zeroth order approximation for the scattering 




(108) 



ift-rrPi = [hi, pi] +Tr 2 {[VL2,/3i2]}. 



(109) 



pw~p w = 
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states in p i2 corresponds to the first order of the Born approximation in the equation 
for pi . 

As we are dealing with bosons we have supplemented the Hartree approximation by a 
bosonic symmetrization procedure, involving the permutation operator P 12 . Note that 
the symmetrization as it was written works only for particles 1 and 2 in orthogonal states: 

l^W) = '">'«> (in) 

as the factor y/2 is the correct normalization factor only in this case. This is almost true 
for a non-degenerate Bose gas. This restriction forces us to treat separately the case in 
which a condensate is present ( T < T c ) . 

We now insert the Hartree-Fock ansatz for p 12 in the hierarchy 3 

»fi|pi = [hup,] + Tr 2 {[lWf/]}. (112) 

In the commutator with Vyi we will encounter 

S(f[ - r2)(l + F 12 ) = (1 + P 12 )6(f[ - f 2 ) = 26(f[ - r 2 ). (113) 

The fact that P i2 commutes with Vu is due to the parity of the delta distribution, 
and P12 acting on a state with two particles at the same position can be replaced by 
the identity. As a consequence, with our zero-range interaction potential, the Fock term 
simply doubles the Hartree term. We finally obtain 

(114) 



E- + u(f) + v(f), Pl 



where V(r ) is the mean field potential 

V(f) = 2g(f\p 1 \r) = 2gp(f). (115) 
The Hartree-Fock Hamiltonian is then 



2 



h HF (l) = ^ + U(r) + 2gp(r). (116) 



3 Note that for the present calculation the regularization of the pseudo-potential is not necessary. 
Indeed by considering plane waves as scattering states in p i2 we suppress any problem of divergences 
in the commutator with V\i , and we can then take V\i as a simple delta distribution. 
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The problem is then formally reduced to the one of an ideal Bose gas moving in a 
self-consistent potential. For g > the mean field corresponds to repulsive interactions, 
as 2gp(f) expels the atoms from the region of high density, while for g < the mean 
field corresponds to attractive interactions. 

4.2.2 Effect of interactions on T c 

Let us now consider the Hartree-Fock one-body density matrix at thermal equilibrium; 
we use the same formula as the ideal Bose gas Eq.(29), replacing hi by the Hartree-Fock 
Hamiltonian: 

h = {exp \P {h HF {l) - fi)} - l}" 1 . (117) 

For k B T AE where where AE is the level spacing of h HF we can perform the 
semiclassical approximation. We obtain for the spatial density as in Eq.(55): 

Psc{r) = T^3/2(zexp[-/?(£/(f ) + 2gp sc (r))}) (118) 

A dB 

At T = T c the argument of # 3 / 2 goes to 1 in the point r m i n where the potential is 
minimal, so that Einstein's condition still holds in the Hartree-Fock approximation: 

ft C (4in)4 = C(3/2). (119) 

For the harmonic trap U(f) = muj 2 r 2 /2 the minimum occurs at the center of the trap, 
^min = so that the chemical potential at the phase transition is given by 

fi = 2gp sc (0). (120) 

It is shifted by the mean field effect with respect to the ideal Bose gas. Using as a small 
parameter p sc (3)g/k B T® , one can derive at constant N [22] the first order change in 
the critical temperature with respect to T c ° , the transition temperature of the ideal Bose 
gas: 

| = -2.5^(0>=-1.33 w ^^/ 6 (121) 

For N = 10 7 atoms of 23 Na in a trap of harmonic frequency u = 2n x 100 Hz, with a 
scattering length a = 50a Bo hr we find T c ° ~ lfi K, and 5T C /T® ~ -2.5 x 10~ 2 , an effect 
for the moment smaller than the experimental accuracy. The fact that ST C is negative 
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for effective repulsive interactions ( a > ) is intuitive: for fixed values of TV and T the 
interacting gas has a lower density at the center of the trap than the ideal Bose gas, so 
that one needs to further cool the gas to get Bose-Einstein condensation. 

• A calculation of ST C beyond mean Held 

The purest situation to study the effect of the interactions on the critical temperature T c 
is the case of atoms trapped in a flat bottom potential; in this case the density is uniform, 
the previously mentioned intuitive mean field effect is suppressed, and our Hartree-Fock 
theory predicts the same critical temperature as the ideal Bose gas. This prediction 
is actually not correct, and rigorous results for the first order correction of T c in ap 1 / 3 
have been obtained recently, by a combination of perturbative theory and Quantum Monte 
Carlo calculations [23]: 



Recent calculations in the many body Green's function formalism confirm this result [24]. 
This effect, if heuristically extended to the trap, is of opposite sign and of the same order 
of magnitude as the mean-field prediction. 

4.3 Hartree-Fock approximation in presence of a condensate 
4.3.1 Improved Hartree-Fock Ansatz 

As already emphasized in the previous subsection the symmetrization procedure of the 
Hartree-Fock prescription Eq.(llO) has to be modified in presence of a condensate. To 
this end we split the one-body density matrix as 



where is the condensate wavefunction, (No) is the mean number of particles in the 
condensate and p[ is the one-body density matrix of the non-condensed fraction. The 
Hartree approximation for the two-body density matrix now reads: 



ST? 



= (2.2 ± 0.25)ap 1/3 + o{ap 1/3 ). 



(122) 



p l = {N )\4>){4>\ + p' 1 



(123) 



Pi ® Pi — (No) 2 |^j <f>\ + remaining Hartree part. 



(124) 
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The first term in the right hand size is already symmetrized; the second term can be 
symmetrized as in Eq.(llO) as it does not involve coexistence of two atoms in the (only) 
macroscopically populated state <j> . We therefore put forward the following Hartree-Fock 
ansatz: 



pi 2 F = (N o ) 2 |0, (/>)((/>, <t>\ + [ ~~~~t=-^ ) remaining Hartree part (-""'" ^ 12 



^_ ^_ j. 

Eliminating the remaining Hartree part with the help of Eq.(124), we finally obtain 



(125) 



~hf (l + Pn \ . ^ (I + P12 
P\i = )Pi ® Pi 



(126) 



y/2 ^ 

In this way we have avoided the double counting of the condensate contribution that 
would have resulted from the prescription Eq.(llO). 



4.3.2 Mean field seen by the condensate 

We replace pu in the first equation of the hierarchy by the improved Hartree-Fock ansatz. 
The first bit of the ansatz gives the same result as in the case T > T c , the second bit 
involves the term: 



TV 2 ([S(fi - r 2 ), I^X^I]) = [|<Kr!)| 2 , \ct>)Ct>\}- 
Splitting pi as condensate and non-condensed contribution we arrive at 



lh dt Pl = 



+ 



EL 

2m 

EL 
2m 



+ U(f) + 2gp(f) - {N )g\4>(r)\ 2 , (N )\4>}(4>\ 



+ U(r)+2gp{r),p' 1 



(127) 

(128) 
(129) 



The non-condensed particles still move in the mean field potential 2gp(f) . On the 
contrary the atoms in the condensate see a different mean field potential: 



2gp(r) - g{N )W)? = ^9P'(r) + s<JV ) IW)| S 



(130) 



where p' is the non-condensed density and (N )\(j)\ 2 is the condensate density. 4 This 

result can be interpreted as follows: An atom in the condensate interacts with non- 

4 A careful reader may argue that we forget here the condition of orthogonality of the eigenstates of 
p[ to (f> . Inclusion of this condition is beyond accuracy of the Hartree-Fock approximation. It will be 
carefully included in the more precise number conserving Bogoliubov approach of §7. 
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condensed particles with the effective coupling constant 2g , and it interacts with another 
particle of the condensate with the effective coupling constant g . 

For repulsive effective interactions ( g > ) this is at the basis of Nozieres' argument 
against fragmentation of the condensate in several orthogonal states: in a box of size L 
in the thermodynamical limit, transferring a finite fraction of condensate particles from 
the plane wave p = to an excited plane wave p = 0(h/L) costs a negligible amount 
of kinetic energy per particle but a finite amount of interaction energy. The transferred 
fraction would indeed be repelled with a stronger amplitude ( 2g rather than g ) by the 
atoms remaining in the condensate. 

4.3.3 At thermal equilibrium 

At thermal equilibrium the one-body density matrix of non-condensed atoms is given 
by the usual Bose distribution for the ideal Bose gas, with the trapping potential being 
supplemented by the mean-field potential: 

Pi = r-r^ 1 n ( 131 ) 

exp {/? + U(f) + 2gp(f) -/*]}- 1 

The condensate wave function has to be a steady state of the total, mean field plus 
trapping potential seen by an atom in the condensate: 

X4>{f) = ~^n^ + [U(f) + 9 {N ° mr )|2 + 2 9p'(r)W)- (132) 

The Hartree-Fock single particle energy A should not be confused with the energy per 
particle in the condensate, as it will become clear in the next section. The occupation 
number of the condensate is related to A by the Bose formula: 

<"°> = ^iri- ( 133 ) 

We now have to solve in a self consistent way the three equations Eq.(131,132, 133). 
In practice, when (N ) is already large, one can assume A = \i , which eliminates one 
unknown A and one equation Eq.(133). 
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4.4 Comparison of Hartree-Fock to exact results 
4.4.1 Quantum Monte Carlo calculations 

The Quantum Monte Carlo method developed by David Ceperley and others allows to 
sample in an exact way the N -body distribution function of a gas of N interacting 
bosons at thermal equilibrium. I.e. the algorithm generates random positions fl,...,r~N 
for the N particles with a probability distribution given by the exact N -body distribu- 
tion function of the atoms. 

On the figure 8 the Hartree-Fock prediction for the radial density of particles in a 
spherical harmonic trap, r 2 p(r) , is compared to the Quantum Monte Carlo result for 
several temperatures below T c . The Hartree-Fock prediction is in good agreement with 
the exact result, except close to T c where it tends to underestimate the number of 
particles in the condensate [25] . 



N(r)/Np 
0.12 - 



f V 



\ o.lo- .} ; 



\0.09|-!.' 



V J 



1.5 2.0 2.5 3.0 3.5 
\ (a) 



N(r) 

0.5 
0.4 
0.3 
0.2 
0.1 
0.0 



0.0 5.0 10.0 15.0 r [a ] 



1 




(b)" 




HF 




/ X 


QMC 




/ 

./ x 






/ 

_ ../ 
1 






0.0 2.0 


4.0 6.0 8.0 


10.0 1 



Figure 8: Radial density of particles, r 2 p(r) , for an interacting Bose gas at thermal equilib- 
rium in an isotropic harmonic trap. Noisy lines: results of a Quantum Monte Carlo simula- 
tion. Smooth solid lines: Hartree-Fock prediction. The curves corresponds to the temperatures 
T/T c ° = 0.88 (a), T/T c ° = 0.7 (b). The number of particles is N = 10 4 and the parameters 
are the ones of 87 Rb. These figures are taken from [25]. 
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4.4.2 Experimental results for the energy of the gas 

At JILA the sum of kinetic and interaction energy of the atoms was measured as function 
of temperature, as we have already explained in §2.3. Whereas the ideal Bose gas model 
was clearly getting wrong for T <T C , the Hartree-Fock prediction [26] is consistent with 
the experimental results over the whole considered temperature range (see figure 9). 




0.5 1 1.5 



T/T 

Figure 9: Expansion energy of the gas per particle and in units of &bT c ° as function of the 
temperature in units of T c ° . The filled rhombi correspond to the experimental results of [11]. 
The straight solid line is the prediction of Boltzmann statistics. The dotted curve is the ideal 
Bose gas prediction. The circles are the numerical solution to the Hartree-Fock approach. The 
curved solid line and the dashed line are approximate solutions to the Hartree-Fock equations. 
The inset is a magnification showing the change of slope of the energy as function of T close 
to T = T c ° . The figure is taken from [26]. 

At very low temperatures ( T < T c /2 ), measurements at MIT have shown that the 
same energy becomes mainly a function of the number of particles N in the condensate. 
By setting p' — in the Hartree-Fock approximation, and using approximations pre- 
sented in the coming section §5.3, an analytical expression can be obtained for the energy, 
in excellent agreement with the experimental results (see figure 5): the energy per particle 
has a power law dependence with N , with an exponent 2/5 , to be contrasted with the 
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constant ideal Bose gas result, and has typical values an order of magnitude larger than 
the zero-point energy of the harmonic oscillator. 



5 Properties of the condensate wavefunction 

In this section we consider the regime of an almost pure condensate, where the non- 
condensed cloud has a negligible effect on the condensate. At thermal equilibrium with 
temperature T this regime corresponds to the limit T <C T c . As we shall see most of the 
experimental results obtained with almost pure condensates can be well reproduced by a 
single equation for the condensate wavefunction, the so-called Gross-Pit aevskii equation, 
derived independently by Gross [27] and Pitaevskii [28]. 



5.1 The Gross-Pitaevskii equation 
5.1.1 From Hartree-Fock 

Let us assume that the density of non-condensed particles is much smaller than the density 
of condensate particles over the spatial width of the condensate: 



(134) 



where N is the mean number of particles in the condensate, <j> is the condensate wave- 
function normalized to unity: 



J d 3 r(/>(r,t)(t>*(r,t) = 1. 



(135) 



In the Hartree-Fock expression of the mean field potential seen by the condensate, derived 
in the previous section §4, we can drop the contribution of the non-condensed particles, 
to get for the evolution of the condensate contribution to the 1-body density matrix: 



ihj t {N,\m\) = 



21 

2m 



+ U(f,t) + gN \<f>(r,t)\ 2 ,N \cl>)(cl>\ 



(136) 



This equation leads to N = constant and to the evolution equation for the condensate 
wavefunction: 



ihd t <f>(r, t) = 



- ^ a + u(f, t) + N gw, t)\ 2 -m 



4>{r,t). 



(137) 
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This non-linear Schrodinger equation is the so-called time dependent Gross- Pit aevskii 
equation. This equation is determined from our Hartree-Fock approach up to an arbitrary 
real function of time, £(t) , as Eq.(136) involves a commutator to which £(t) does not 
contribute. In general the precise value of is considered as a matter of convenience, 
as it can be absorbed in a redefinition of the global phase of <j> . The knowledge of 
the value of £(t) can become important when one is interested in the evolution of the 
relative phase of two Bose-Einstein condensates. The value of has been derived in 
[29] assuming a well defined number of particles in the condensate. If the condensate is 
assumed to be in a Glauber coherent state that is a quasi-classical state of the atomic 
field with a well defined relative phase (see §8) one obtains £(t) = as we will see in 



When the gas is at thermal equilibrium, the only time dependence left for is a 
global phase dependence. The most convenient choice is to assume d t (f> = so that £(t) 
is a constant. As shown in §4.3.3 this constant is very close to the chemical potential 
of the gas as N Q is large so that we get the so-called time independent Gross- Pit aevskii 
equation: 



Both the time independent and the time dependent Gross-Pitaevskii equations can be 
solved numerically. But, as explained in the next part of this section, the fact that the 
trap is harmonic allows one to find very good approximate analytical solutions. 

5.1.2 Variational formulation 

Variational calculus turns out to be a very fruitful approximate technique in the solution 
of the Gross-Pitaevskii equation. We therefore derive here a variational formulation of 
the Gross-Pitaevskii equation. 

• Time independent case 

The time independent Gross-Pitaevskii equation can be obtained from extremalization 
over <j> of the so-called Gross-Pitaevskii energy functional: 



5.1.3. 



MO= -^A + ^(f) + V <#(f)| 2 



(138) 



E[(j>,(f>*}=N J Sr 



2m 



gmd0| 2 + U{r)W)\ 2 + \N,gW)\ A (139) 
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with the constraint that (/> is normalized to unity. 

Proof: We take into account the normalization constraint with the method of La- 
grange multiplier, so that we simply have to express the fact that <j> extremalizes without 
constraint the functional: 



X[4>,<t>*] = E[<t>,<j>'} - XN J d 3 f 4>(fW(r)- 



(140) 



The parameter A is the Lagrange multiplier. We calculate the first order variation of X 
due to an infinitesimal arbitrary variation of the condensate wavefunction: 



We obtain: 
SX = N J d 3 r 



2m 



(j>(r) -> 0(f) + <ty(f). (141) 



grad^*-grad0 + ^(f)^V + A r o^V> 2 -A^> + c.c. (142) 



We modify the variation of the kinetic energy term by integrating by part, assuming that 
(f> vanishes at infinity: 

J d 3 f (grad<ty* • grad0 + c.c.) = - J d 3 f (<ty*A0 + c.c). (143) 

The variation SX has to vanish for any 5$ . We can take as independent variables the 
real part and the imaginary part of 8(j) , or equivalently 8(j) and 5(f)* as it amounts 
to considering independent linear superpositions of the real and imaginary part. We 
therefore obtain: 



-^A + t/(f) + ^|0(f)| 2 -A 



(p(f) = 0. 



(144) 



We recover the time independent Gross- Pit aevskii equation, with A = \i , which gives a 
physical interpretation to the Lagrange multiplier A . 

• Time dependent case 

The time dependent Gross- Pit aevskii equation with the choice £(t) = is obtained over 
a time interval [^1,^2] from extremalization of the action: 

rt2 



A 



= / 



dt 



f (^|||0>-c.c.)jVo-W),^(*)] 



(145) 



with fixed values of \<j>(t = ti)) and \<j>(t = t 2 )) . 
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• Physical interpretation of the Gross-Pitaevskii energy functional 

We now show that E[</>, </>*] is simply the mean energy of the gas in the Hartree-Fock 
approximation in the limit of a pure condensate. As the N -body Hamiltonian is a sum 
of one-body and two-body (binary interaction) terms, 

H = j:hi + WVi3 (146) 
i=i z i^j 

the mean energy of the gas involves the one-body and two-body density matrices: 

(H) = Ti[hipi] + ^Tr[y 12 p 12 ]. (147) 
In the limit of a pure condensate we keep only the condensate contribution to pi : 

Pi-iVol^l (148) 
and we approximate p i2 by the Hartree ansatz 

pi 2 ~pi<g>pi. (149) 

We then obtain E[(j),(j)*] = (H) . It was actually clear from the start that E[$,<f>*] was 
the sum of kinetic energy, trapping potential energy and mean field interaction energy of 
the condensate. 

A different and interesting point of view at zero temperature is to use directly a 
Hartree-Fock ansatz for the ground state wavefunction of the gas, assuming that all 
the particles are in the condensate: 

|*) = |iV:0> = |0><8)...<8)|^>. (150) 

The mean energy of |\I>) for the interaction potential g6(fi - r 2 ) is then 



r *2 



(151) 



which differs from Eq.(139) in the limit N = N only by the occurrence of a factor (N—l) 
rather than N in front of the coupling constant g , ensuring that the interaction term 
disappears for N = 1 ! 
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• What is the chemical potential ? 

At zero temperature, assuming a pure condensate N ~ N , the usual thermodynamical 
definition of the chemical potential \i reduces to: 

»=^m*w^ (152) 

where we have made appear the explicit dependence of E on N Q . When one takes the 
total derivative of E with respect to N Q , one gets in principle a contribution from the 
implicit dependence of E on N through the N Q dependence of <f>, (j>* ; actually this 
contribution vanishes as the variation of E due to a change in </5>, <f>* vanishes to first 
order in this change. We therefore get 



dN ' UJ 8N 

r fc2 



= I d 3 r 



(153) 



This quantity coincides with the chemical potential indeed, as can be checked by multi- 
plying the time independent Gross-Pitaevskii equation by <f>* and integrating over the 
whole space. As g does not have the factor 1/2 in Eq.(153), whereas it is multiplied by 
1/2 in the expression for E[(j), cjf] , we see that in the interacting case g ^ : 

(154) 

that is the chemical potential \i differs from the mean energy per particle. 

5.1.3 The fastest trick to recover the Gross-Pitaevskii equation 

Starting from the second quantized form of the Hamiltonian, 

H = jdX ^(l)M(l) + \ j dlf d2 ^(1)^(2)^(2)^(1) (155) 

where 1 and 2 stand for three-dimensional coordinates in real space, one first derives 
the Heisenberg equation of motion for the field operator: 

ihj t m = im,H] = d^ Hl) H (lse) 

= M(l) + |d2^ t (2)y 13 ^(2)^(l) (157) 
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and then replaces the quantum field operator by a classical field: 




As V12 is the pseudo-potential, the equation that we get for (j> is the time dependent 
Gross-Pitaevskii equation with f (£) = . 

This sheds a new light on the Gross-Pitaevskii equation: the Gross-Pitaevskii equation 
is the equation of motion of the atomic field in the classical approximation, neglecting 
quantum fluctuations of the field. A Bose-Einstein condensate is a classical state of the 
atomic field, in a way similar to the laser being a classical state of the electromagnetic 
field. 



5.2 Gaussian Ansatz 

In this subsection we look for a variational solution to the Gross-Pitaevskii equation in a 
harmonic trap, using a Gaussian ansatz for [30]. The choice of a Gaussian is natural 
in the non-interacting limit g — > , where it becomes exact. It turns out to give also 
interesting results in presence of strong interactions. 



5.2.1 Time independent case 

Consider for simplicity an isotropic harmonic trap, where the atoms have the oscillation 
frequency uj . We assume the following Gaussian for the condensate wavefunction: 

the spatial width o being the only variational parameter. The mean energy per particle 
can be calculated exactly for this ansatz: 

_ Efap] 3h 2 3 2 2 h 2 N Q a 1 

e = — — = + -muj a H — (160) 

AT 4ma 2 4 ma 3 J2n 



The form of the result is intuitive: the kinetic energy term scales as Ap 2 x , where Ap x = 
h/(2Ax) — h/(\/2a) ; the trapping potential energy scales as a 2 and the interaction en- 
ergy per particle is proportional to the coupling constant g — A7rh 2 a/m and to the typical 
density of atoms in the gas, N Q /a 3 . Taking the harmonic oscillator length (h/muj) 1 / 2 
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as a unit of length and the harmonic quantum of vibration Huj as a unit of energy we 
get the simple form: 

" + h (161) 



3 

e = - 



4 [a 2 

where the only physical parameter left is 



1 2 



*=i/H!r=. (162) 

Jn/muj 

This parameter x measures the effect of the interactions on the condensate density: The 
case x < 1 corresponds to the weakly interacting regime, close to the ideal Bose gas 
limit x — ; the case x ^> 1 corresponds to the strongly interacting regime. 

• case a > 

In the case of effective repulsive interactions between the particles, the dependence of e 
with a is plotted in figure 10. In the limit a — > , the energy e is dominated by the 
positively diverging repulsive interaction ( ~ 1/<t 3 ). For large a the trapping potential 
term ~ a 2 dominates. The function e has a single minimum, in a = <j , solving 

^(a ) = 0^a 5 = a + x- (163) 

For x <C 1 one recovers the ground state of the harmonic trap, with a Q — 1 . For x ^ 1 
the condensate cloud becomes much broader than the ground state of the harmonic trap, 

a Q ~ x 1/5 (x 7V 1/5 . (164) 

In this regime one can check that the kinetic energy term becomes negligible as compared 
to the trapping energy: 

Eki » 1 ~ — (165) 



Strap <? 4 X 4/5 

so that the steady state of the condensate is an equilibrium between the trapping potential 
and the repulsive interactions between particles. This regime will be studied in detail in 
the next subsection. 



• case a < 
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Figure 10: Energy per particle in the condensate in units of Tiuj as function of the variational 
width g in units of (h/mu) 1 / 2 . Case of effective repulsive interactions a>0. 

For effective attractive interactions between the particles the shape of e as function of 
a depends on the balance between kinetic and interaction energy (see figure 11). The 
interaction energy is negatively diverging as o — > always faster than the positively 
diverging kinetic energy so that a — is always a minimum of e , with e — — oo : 
the condensate is in a spatially collapsed state ! Of course the Gross-Pitaevskii equation 
not longer applies for a too small a , as the validity of the Born approximation requires 
k\a\ ~ |a|/a^C 1 ■ For \x\ larger than some critical value \xc\ , this collapsed minimum 
is the only one of e(a) so that we do not find any stable solution for the condensate 
wavefunction. When \x\ is smaller than |% c | the kinetic energy term, which is opposed 
to spatial compression of the gas, is able to beat the attractive energy over some range of 
a , so that a local minimum of s(a) appears, in a = a Q , separated from the collapsed 
minimum by a barrier. 

To calculate |x c | we express the fact that the stationary point of e in a = <To has 
now a vanishing curvature (inflexion point of e ) : 

\ / O"=0"o 
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Figure 11: For an effective attractive interaction a < between the particles energy per 
particle in the condensate in units of hu as function of the variational width a in units of 
(h/mou) 1 / 2 . The curve has two possible shapes, (a) with two minima when \x\ is smaller than 
a critical value \xc\ , an d (b) with a single minimum for |x| > |Xc| ■ 
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X=Xc 



= 



\da 2 

\ / <J = <JQ 

By eliminating a between these two equations we obtain 

4 



Xc = 



"5 5 / 4 



= -0.5350... 



(167) 



(168) 



This result can be rephrased in terms of a maximal number of atoms iVJ that can be 
put in the condensate without inducing a collapse, according to a Gaussian ansatz: 

|Xe| - -0.67. (169) 



N$\a\ 



A more precise result has been obtained by a numerical solution of Gross- Pit aevskii 
equation, not restricting to the subspace of Gaussian wavefunctions [31]: no solution of 
the time independent Gross-Pit aevskii equation is obtained for N > Nq , where 



~ -0.57. 



(170) 



muj 
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By a generalization of the Gaussian ansatz to the case of a non-isotropic harmonic 
trap one can also get a prediction of Nq for the parameters of the lithium experiment 
of Hulet's group [32]. In the experiment the traps frequencies are uj z = 2tt x 117 Hz and 
&x,y — 27T x 163 Hz, and the scattering length is a — -27 Bohr radii. The Gaussian 
prediction is then iVJ ~ 1500 , consistent with the experimental results. 

• Physical origin of the stabilization for a < 

In a harmonic trap, the energy of the ground state level is separated from the energy 
of excited states by Huj . At low values of x the mean interaction energy per particle, 
~ p\g\ , where p is the density, is much smaller than Huj so that it cannot efficiently 
induce a transition from the ground harmonic level to excited harmonic levels. Initiation 
of collapse on the contrary requires that the wavefunction <j> expands on many excited 
levels in the trap, so that the density \(j)\ 2 can exhibit a high density peak narrower than 
muj . We therefore intuitively reformulate the non-collapse condition as 

p\g\ < fiuj. (171) 

Estimating p as Nq/ '(/i/mcj) 3 / 2 we recover a Nq scaling as ijh/muj/\a\ . This reasoning 
also applies to the gas confined in a cubic box with periodic boundary conditions, as we 
shall see in section §6 of the lecture. 




5.2.2 Time dependent case 



As done in [33, 34] the Gaussian ansatz can be generalized to the time dependent case. 
We assume here for simplicity that the condensate, initially in steady state, is excited only 
by a temporal variation of the trap frequencies u a (t) ; then no oscillation of the center 
of mass motion of the condensate takes place, <j> remaining of vanishing mean position 
and momentum. The Gaussian ansatz then contains only exponential of terms quadratic 
with position, its does not involve exponential of terms linear with position: 

e iS(t) 



^ 4 [rw*)l 



1/2 



exp 



-E 



(172) 



We do not assume that the trap is isotropic, so we have as variational parameters 3 spatial 
widths a a (a = x,y,z), 3 factors 7 a governing the spatially quadratic phase and a 
global phase S . 
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One gets time evolution equations for the variational parameters by inserting the 
ansatz for in the action A of Eq.(145) and by writing the Lagrange equations ex- 
pressing the stationarity condition. It turns out that j a can be expressed in terms of 
the widths and their time derivatives: 

7q = "as" (173) 

so that one is left with equations for the a a 's. Taking u~ l as a unit of time, \jfijmuj 
as a unit of length, where uj is an arbitrary reference frequency, we get: 

0-o+^ff„ = 4 + ( 174 ) 

where the trap frequencies are uj a = v a uo and x is defined in Eq.(162). In the absence 
of interaction ( x — ) these evolution equations become exact, and give a remarkable 
(and known !) result for the time dependent harmonic oscillator. In the interacting case 
( X ^ ) these equations can be cast in Hamiltonian form as the "force" seen by the 
variable a a derives from a potential. The corresponding dynamics is non linear and non 
trivial; chaotic behavior has been obtained in [35] in the limiting regime of x ^ 1 where 
the l/<7^ can be neglected. 

One can use Eq.(174) to study the response of the condensate to a weak excitation, 
the trap frequency Lo a in the experiments being typically slightly perturbed from its 
steady state value cj a (0) for a finite excitation time. Linearizing the evolution equations 
in terms of the deviations of the a 's from their steady state value: 

a a (t) = o% + 5a a (t) (175) 

one gets a three by three system of second order differential equations for the 6 a 's. 
Looking for eigenmodes of this system, one finds three eigenfrequencies [34]. Their values 
have been compared to experimental results at JILA [36], see Fig. 12: the agreement is 
very good, not only in the weakly interacting regime x ^ 1 but also in the regime 
X ^> 1 , where the Gaussian ansatz for the condensate wavefunction has no reason to be 
a good one! The explanation of this mystery is given in §5.4.1. 

5.3 Strongly interacting regime: Thomas-Fermi approximation 

In this subsection we focus on the strongly interacting regime: the scattering length is 
positive, with the dimensionless parameter x °f Eq.(162) much larger than one. This 



The condensate wavefunction 



55 




2 



o 



o 



o 



1.4 - 



1.2 







2 



4 



6 



10Xv r ' 



1/2 



Figure 12: Frequencies of two eigenmodes of a condensate in a cylindrically symmetric harmonic 
trap, in units of the radial trap frequency v r , as a function of a parameter proportional to \ 
measuring the strength of the interactions. Plotting symbols: measurements at JILA [36]. Solid 
lines : predictions of the Gaussian Ansatz [34]. 

regime is the so-called Thomas- Fermi regime. As we now see analytical results can be 
obtained in this limit. 

5.3.1 Time independent case 

If we put a large enough number of particles into the condensate the atoms will experience 
repulsive interactions that will increase the spatial radius of the condensate to a value R 
much larger than the one of the ground state of the harmonic trap: 



For increasing value of A^ , R increases so that the momentum width of the condensate, 
scaling as Ti/R as 4> Q is a non-oscillating function of the position, is getting smaller and 
smaller. More precisely we find that the typical kinetic energy of the condensate becomes 




(176) 
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much smaller than the typical harmonic potential energy of the condensate: 

jC/kin mi? 2 \ '1 \ ^ 1 (177) 




The mechanical equilibrium of the condensate in the trap then comes mainly from the 
balance between the expelling effect of the repulsive interactions and the confining effect 
of the trap. 

In this large R regime we neglect the kinetic energy term in the Gross-Pitaevskii en- 
ergy functional, which leads to a functional of the condensate density only (similarly to the 
Thomas- Fermi approximation for electrons). This approximation amounts to neglecting 
the A(f> term in the Gross-Pitaevskii equation: 

ii^r) ~ U(f )(j>{r ) + N*gW ). (178) 

Taking <f> to be real we find that 

in the points of space where fi > U(f) , otherwise we have 4>(r) = . 

This very important, Thomas-Fermi result Eq.(179) can also be obtained in a local 
density approximation point of view. A spatially uniform condensate with a chemical 
potential \i and in presence of a uniform external potential U has a density iV 1 <^ | 2 = 
(li — U)/g. Applying this formula with a r dependent potential U gives again Eq.(179). 
A local density approximation can be used only if the density of the condensate varies 
slowly at the scale of the so-called "healing length" £ , introduced in §5.3.4; one can check 
that the condition £ < R is indeed satisfied in the Thomas-Fermi regime. 

We specialize Eq.(179) to the case of a harmonic but not necessarily isotropic trap: 

U{r) = \m^lrl (180) 

where a — x,y,z label the eigenaxis of the trap. The boundary of the condensate 
\i— U(f) is then an ellipsoid with a radius R a along axis a given by: 

M=imwX- ( 181 ) 
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The condensate wavefunction can be rewritten in terms of these radii: 

= (ji_\ 1/2/ ~ 2Nl/2 



N gJ 



('-?!) 



(182) 



Using the normalization condition of <j> to unity we can also express the "normalization" 
factor ^fi/N Q g in terms of the radii. The integral of |0| 2 can be calculated in spherical 
coordinates after having made the change of variable u a = r a /R a . This leads to 

/ \ 1/2 / \ 

Ml 



,N g) 



15 



8tt l[R a 



(183) 



Eliminating R a in terms of \i thanks to Eq.(181) we can calculate the chemical potential: 

1 2/5 



15- 



N a 



(h/mti)) 1 / 2 _ 

where Co is the geometrical mean of the trap frequencies: 

L) = ((Jj x (jjy(jj z ) 1/3 . 

We can now see that in the limit % > 1 the chemical potential fi satisfies 



(184) 



(185) 



(186) 



which is a convenient way of defining the Thomas- Fermi regime. 

We can now compare these Thomas-Fermi predictions to the MIT experimental results 
on the energy of the condensate [12]. In the experiment the trapping potential is switched 
off abruptly, so that the energy of the gas abruptly reduces to £^ red = E kin + £^ int ~ E- mt ; 
afterwards the cloud ballistically expands, E- mt is converted in kinetic expansion energy 
that can be measured. In the Thomas-Fermi approximation the integral of 7Vq #|</5>| 4 /2 
can be done, which leads to 

£mt ^ fan oc AT 7/5 . (187) 

The resulting dependence in A^ is in good agreement with the MIT results, see Fig. 5. 

From the expression of the chemical potential we can also calculate the total energy 
of the condensate in the trap , as \i = d No E : integrating over gives 



E ~ -fiN . 

One can then check explicitly that fi ^ E/N \ 



(188) 
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5.3.2 How to extend the Thomas- Fermi approximation to the time dependent 
case ? 

We would like to analyze time dependent situations encountered in the experiments, e.g. 

• ballistic expansion of the gas: this is a crucial example, as it is a standard experi- 
mental imaging technique of the condensate 

• collective excitations: response of the condensate to a modulation of the trap fre- 



in the strongly interacting regime. An immediate generalization of the Thomas-Fermi 
approximation consisting in neglecting the kinetic energy of the condensate is now too 
naive! In the case of ballistic expansion for example the interaction energy is gradu- 
ally transformed into kinetic energy when the cloud expands so kinetic energy becomes 
important! 

The trick is actually to split the kinetic energy in two contributions, one of them 
remaining small and negligible in the time dependent case. This is performed using the 
so-called hydrodynamic representation of the condensate classical field, split in a modulus 
and a phase: 



where S has the dimension of an action and p is simply the condensate density. The 
mean kinetic energy of the condensate then writes 



As we shall see during ballistic expansion of the condensate the density p remains a 
smooth, slowly varying function of the position so that it has a very small contribution to 
the kinetic energy; most of the kinetic energy induced from interaction energy is stored 
in the spatial variation of the phase of the condensate wavefunction. 



quencies 



(189) 




(190) 
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5.3.3 Hydro dynamic equations 



In this subsection we rewrite the time dependent Gross-Pitaevskii equation in terms of 
the density p and the phase S . This can be done of course by a direct insertion of 
Eq.(189) in the Gross-Pitaevskii equation. 

A more elegant way is to use the covariant nature of the Lagrangian formulation of 
the Gross-Pitaevskii equation, Eq.(145). We rewrite the density of Lagrangian in terms 
of p and S : 



£= - 



pdtS + £ (^/p) 2 + p^p- + u <?> $p + § p 2 



(191) 



(192) 



An evolution equation for an arbitrary coordinate Q(f,t) of the field is obtained from 
the Lagrange equation: 

/ ac \ „ / ac \ = dc 

' \d(d t Q) ) + ^ r ° \d(d Ta Q) ) dQ- 

We first specialize the Lagrange equations to the choice Q — ^fp ; dividing the result- 
ing equation by 2-^/p we obtain 

dt S + ^l AS )^U + P9 =^f. (193) 
Then we set Q — S in the Lagrange equations which leads to 



d t p + div 



^gradS 
Im 



= 0. 



(194) 



This last equation looks like a continuity equation. This is confirmed by the following 
physical interpretation of grad S . It is known in basic quantum mechanics that the 
probability current density associated to a single particle wavefunction <f> is 



Jproba 



= ^[^ ad< ^- cx -' ' 



2im 



(195) 



Multiplying this expression by , as there are particles in the condensate, and 
introducing the (p, S) representation of (j) we get the following expression for the current 
density of condensate particles: 



grad S 
J = P = pv 



m 



(196) 
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where v is the so-called local velocity field in the gas. 

Equation (194) is therefore the usual continuity equation: 

d t p + div[pv] = 0. (197) 

The other equation (193) can be turned into an evolution equation for the velocity field 
by taking its spatial gradient: 

t2 



md t v + grad 



= 0. (198) 



This looks like the Navier-Stockes equation used in classical hydrodynamics, in the 
limiting case of a fluid with no viscosity. The term grad(|rra; 2 ) looks unusual but using 
the fact that v is the gradient of a function S/m one can put it in the usual form of a 
convective term: ^ 

grad (^-mv 2 ^ — m(v • grad)i7. (199) 

A difference with classical hydrodynamics is the so-called quantum pressure term, involv- 
ing 

-- 200 

the only term in the equations (197,198) where h appears. 

5.3.4 Classical hydrodynamic approximation 

The classical hydrodynamic approximation consists precisely in neglecting the quantum 
pressure term Eq.(200) in the equation (198) for the velocity field of the condensate. 

We can estimate simply the validity condition of this approximation. Denoting d 
a typical length scale for the variation of the condensate density p(f) we obtain the 
estimate 

Comparing the quantum pressure term Eq. (200) to the classical mean field term pg yields 
the condition 

ti 2 

< gp(f) - p max £ (202) 



md 2 
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where yO max is the maximal density (usually at the center of the trap). This validity 
condition can be reformulated in terms of the healing length, 



\2mp max gJ 



Note that f is sometimes also called coherence length, which can be confusing. 

Why this name of healing length for f ? Imagine that you cut with an infinite wall a 
condensate in an otherwise uniform potential. Right at the wall the condensate density 
vanishes; far away from the wall the density of the condensate is uniform. The condensate 
density adapts from zero to its constant bulk value over a length typically on the order 
of f . This can be checked by an explicit solution of the Gross-Pitaevskii equation: 



where z = is the plane of the infinite wall. This explicit solution shows that at 
a distance z ^> f from the infinite wall there is no more any effect of the boundary 
condition <j>(x,y,z = 0) = . This is to be contrasted with the case of the ideal Bose 
gas: the ground state between infinite walls separated by the length L then scales as 
sin(7rz/L) and depends dramatically on L . 

For a moderate excitation of the condensate by a modulation of the trap frequencies, 
or in the course of ballistic expansion of the condensate, we shall see that the only typical 
length scale for the variation of the condensate density is the radius R of the condensate 
itself. One can then check that in the Thomas-Fermi regime the classical hydrodynamic 
approximation indeed applies: 



In the Thomas-Fermi regime we therefore neglect the quantum pressure term to obtain 




(204) 




(205) 



m d t + (v • grad) j v(r, t) = -grad (U (r, t) + gp(f, t)) = F(f, t). 



(206) 



This equation is then a purely classical equation, Newton's equation in presence of the 
force field F written in Euler's point of view. The operator between square brackets is 
simply the so-called convective derivative. 
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It is instructive to rewrite Eq.(206) in Lagrange's point of view. One then follows a 
small piece of the fluid in course of its motion. Denoting f(t) the trajectory of the small 
piece of fluid we directly write Newton's equation: 

m % f{t) = = ~ (c/ + 9p) ] *)■ (207) 

This equation automatically implies the continuity equation (197) and the Euler equation 
(206). The unusual feature is that the force field depends itself on the density of the 
gas, so that we are facing here a self-consistent classical problem, corresponding formally 
to the mean field approximation for a collisionless classical gas! A surprising conclusion, 
knowing that we are actually studying the motion of a Bose-Einstein condensate! 



5.4 Recovering time dependent experimental results 
5.4.1 The scaling solution 

It turns out that the self-consistent classical problem Eq.(207) can be solved exactly for 
the particular conditions of a gas initially at rest and in a harmonic trap. 

At time t = we assume a steady state Bose-Einstein condensate in the trap, of 
course in the Thomas-Fermi regime so that the classical hydrodynamic approximation 
is reasonable. The steady state of Eq.(207) corresponds to a force field F vanishing 
everywhere, so that 

U(f) + gp(f ) = constant. (208) 

One recovers the stationary Thomas-Fermi density profile, the constant being determined 
from the normalization condition of p and therefore coinciding with the Thomas-Fermi 
approximation for \i . 

At time t > the trapping potential remains harmonic with the same eigenaxis [37] 
but the eigenfrequencies of the trap can have any time dependence: 

U(r,t) = l £ ™&ty»- (209) 

z a=x,y,z 

Then any small piece of the fluid with initial positions r a (0) along axis a will move 
according to the trajectory 

r a (t) = X a (t)r a (0) (210) 
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where the scaling factors X a (t) depend only on time, not on the initial position of the 
small piece of fluid. In other words the density of the gas will experience a mere (possibly 
anisotropic) dilatation 

^ t)= m^A { ^) ht=0 )- (211) 

We can see simply why the ansatz Eq.(210) solves indeed Eq.(207) for a harmonic trap. 
As the initial density in the trap has a quadratic dependence on position, so will have the 
density at time t . The gradient — gvad(pg) appearing in the expression of the force field 
will then be a linear function of the coordinates; so is the harmonic force — grad U(f, t) . 
Newton's equation is therefore linear in the coordinates; dividing it by r a (0) one then 
gets equations for X a (t) irrespective of the initial coordinates r a (0) ! 

More details are given in [38, 39], we give here the equations for the scaling parameters: 

^rM t )= m°m -"lm a (t), a = x,y,z (212) 

at ^a A x A y A z 

with the initial conditions 

A Q (0) = (213) 
|A Q (0) = (214) 

since the condensate is initially at rest. 

Finally we make the connection between these scaling solutions and the equations for 
the spatial widths a a obtained in Eq.(174) from a time dependent variational Gaussian 
ansatz for the condensate wavefunction. We are here in the Thomas- Fermi regime x ^> 1 
so that the terms can be neglected in Eq.(174). The steady state solutions for the 

a a 's are then <7 a (0) ~ x^^V^O) where v is the geometrical mean of the initial 
frequencies ^ a (0) , and the quantities cr a (t)/a a (0) then obey the same equations as the 
A a 's! The Gaussian ansatz, which has the wrong shape in the Thomas-Fermi regime, is 
however able to capture the right scaling nature of the solution! This explains why the 
collective mode frequencies obtained from Eq.(174) are a good approximation, not only 
when % ^ 1 ? where the Gaussian ansatz was expected to hold, but also in the strongly 
interacting regime x ^> 1 ■ 
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5.4.2 Ballistic expansion of the condensate 

At time t = + the trapping potential is turned off suddenly. The scaling parameters 
then satisfy the simpler equations 

>>^- 

These equations are still difficult to solve analytically. In the experimentally relevant 
regime of cigar-shaped traps, with u z (0) <C uj x (0) = u y (0) , one can find an approximate 
solution [38]. 

Experimentally the scaling predictions have been tested carefully. First one can see if 
the ballistically expanded condensate density has indeed the shape of an inverted parabola 
[38]. Second one can measure the radii of the condensate as function of time to see if 
they fit the scaling predictions [40]. Both tests confirm the scaling predictions in the 
Thomas- Fermi regime. 

5.4.3 Breathing frequencies of the condensate 

A typical excitation sequence of breathing modes of the condensate proceeds as follows. 
One starts with a steady state condensate in the trapping potential. Then one modulates 
one of the trap frequencies for some finite time t exc , at a frequency close to an expected 
resonance of the condensate. Then one lets the excited condensate evolve in the unper- 
turbed trap for some adjustable time t osc . Finally one can perform imaging of the cloud, 
e.g. by performing a ballistic expansion of the condensate and measuring the aspect ratio 
of the expanded cloud. By repeating the whole sequence for different values of t osc one 
can reconstruct the aspect ratio as function of t osc ■ 

Such a procedure has been used at JILA and at MIT. In figure 13 are shown results 
obtained at MIT in a cigar-shaped trap, for a modulation of the trap frequency along 
the slow (that is weakly confining) axis z . The solid line corresponds to the prediction 
of scaling theory, the input parameters being (i) the oscillation frequencies uj a of the 
atoms in the trap, (ii) the precise way the excitation is performed, and (iii) the duration 
of ballistic expansion; the agreement between theory and experiment is good, considering 
the fact that there is no fitting parameter [38]. 

If one is interested only in the frequencies of the breathing eigenmodes of the conden- 
sate it is sufficient to linearize the equations of the scaling parameters around their steady 
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state value: 



dt 2 



5X a = -2uj 2 a (0)5X a -uj 2 a {0)J2^ 



(216) 



and find the eigenvalues of the corresponding three by three linear system (see also §6.3.3). 
For a trap with cylindrical symmetry one gets the eigenfrequencies O = \/2u±(Q) and 



n 2 = ; 



3wf(0) + 4wl(0) ± (9w*(0) + 16wi(0) - 16wf (O)wi(O)) 



1/2 



(217) 



The mode observed at MIT corresponds to the - sign in the above expression; as the 
trap was cigar-shaped in the experiment, uj±(0) ^> (jo z (0) so that one has the approximate 
formula 

, 5 x 1/2 



/5\ i/z 



(218) 



0.35 - 



o 

CD 
Q. 

3 




0.25 



100 



t 0EC [ms] 



Figure 13: Aspect ratio of the excited and ballistically expanded condensate as a function of 
free oscillation time i sc • The expansion time is 40 ms, the unperturbed trap frequencies are 
kU_(0) = 27r x 250 Hz, u z (0) = 2n x 19 Hz. Solid line: theory. Diamonds: experimental data 
obtained at MIT. 
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6 What we learn from a linearization of the Gross- 
Pitaevskii equation 

There are several important motivations to perform a linearization of the Gross-Pitaevskii 
equation around a steady state solution O : 

• as the Gross-Pitaevskii equation is a non-linear equation it is crucial to check the so- 
called "dynamical" stability of the steady state solution. More precisely one has to 
check with a linear stability analysis that any small deviation 5(f> of the condensate 
wavefunction from </5> does not diverge exponentially with time. Otherwise </5> 
may not be physically considered as a steady state as even very small perturbations 
will eventually induce an evolution of the condensate wavefunction far from (f> Q . 

• as a byproduct of linear stability analysis we obtain a linear response theory for the 
condensate very useful to interpret experiments which apply a weak perturbation 
to the condensate. 

• another important byproduct is the Bogoliubov approach which gives a description 
of the state of the non-condensed particles that is still approximate but more ac- 
curate at low temperature (typically k B T < fi) that the Hartree-Fock approach. 
This allows to check the so-called "thermodynamical stability" of the condensate 
and will be the subject of §7. 

6.1 Linear response theory for the condensate wavefunction 

6.1.1 Linearize the Gross-Pitaevskii solution around a steady state solution 

Let 4>o(r) be a steady state solution of the Gross-Pitaevskii solution in the time inde- 
pendent trapping potential Uo(f) : 



The trapping potential is then slightly modified by a time dependent perturbation 
5U (r , t) , resulting in a total trapping potential 



0= -— A + £/o + ^o|0o| 2 -M 0o- 



(219) 



U{?,t) = U Q {r f ) + 6U{?,t). 



(220) 
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The condensate wavefunction, initially equal to </5> , evolve according to 



ihd t (j) = 



-^A + U + gN W 2 -,x 



(221) 



As 8U is so small we assume that (j) experiences only a small deviation from O : 

^(f,t) = ^(f) + <ty(r,f) (222) 

so that we can linearize Eq.(221) in terms of 8<j> . Neglecting the second order product of 
8(j) and 8U we obtain: 



ihd t S(f> = 



-—A + U -ii 
2m 



5ct> + ZgN^lfadt + gNz$84>* + SUfa. (223) 



Note the presence of the factor 2 in front of the term proportional to g8(j) ; it turns out 
(and this should become clear in the Bogoliubov approach) that this factor 2 has the same 
origin as the one in the Hartree-Fock potential Eq.(115) for the non-condensed particles. 
As (j) remains normalized to unity, as 0o was, we note that to first order in S(f> , 



J d 3 f [8^t)^(f) + UrWim = 0. 



(224) 



A peculiar feature of Eq.(223) is that, though it is obtained from a linearization pro- 
cedure, it is not a linear equation for 5(j) in the strict mathematical sense: if 5(j) is a 
particular solution of the homogeneous part of this equation (set SU = ), the function 
aS(f> (where a is a constant complex number) is generally not a solution of the homoge- 
neous part anymore because of the coupling of 8(f) to 8(f)* . There are several possibilities 
to restore this linearity. A first one is to consider as unknown functions the real part 
and the imaginary part of 6<j> . A second, more elegant method, more common in the 
literature, is to introduce formally as unknown the two-component column vector: 



^(f, t) 
6</>*(r,t) 



(225) 



the functions 8<j) and 8<j)* being now considered as independent. We then rewrite 
Eq.(223) as the linear system: 



ihdt 



8cf>(r,t) 
5<f>*(r,t) 



GP 



6<f>(f,t) 



+ 



S(f,t) 
-S*(r,t) 



(226) 
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with a source term S(f,t) = 8U(r,t)(j)Q(r) and a linear operator 

r _ ( H GP + gN \4> \ 2 gNofi \ 

Lgp -[ - 9 n 4>? -[h gp + 9 nM}' ) (227) 

where we have introduced the Gross-Pitaevskii Hamiltonian: 

h 2 

H GP = -—A + U Q + gNolfal 2 - fi. (228) 
2m 

Note the presence of complex conjugation in the second line of Cqp ; it also applies to 
the potential U , without effect here as U , hermitian function of r , is real; it should 
not be forgotten if situations where the potential contains a complex term such as — Q,L Z 
where L z is the angular momentum operator (inertial term in a frame rotating at angular 
velocity Q ). 

As the operator Cqp is time independent the general method to determine the time 
evolution of 5(f) is to diagonalize Cqp an d expand 8(j) on the corresponding eigenmodes. 
At this stage one faces a slight difficulty: it turns out that C GP is not diagonalizable, 
that is the set of all eigenvectors of Cqp does not form a basis (in general one vector 
is missing to span the whole functional space). Mathematically this can be solved by 
putting C GP into the so-called Jordan normal form. Here we use the more physical 
following procedure. 

6.1.2 Extracting the "relevant part" from S(j) 

We split 8(j> on a component along O and a part orthogonal to O : 

6(/>(f,t) = v(t)Mr) + <ty±(r,t). (229) 

From Eq.(224) valid to first order in 5(f> we realize that rj = (<^ |^) is such that rj(t) + 
rj*(t) = , so that rj(t) is purely imaginary and can be reinterpreted as a change of phase 
of 0o : 

t) ~ e^Vo(r ) + 5</> ± (?,t). (230) 

One then sees that this change of phase has no consequence on the one-body density 
matrix of the condensate, up to first order in 8(f) : 

|^| ~ |<^o><^o| + |^o> (S4>±\ + \6ti_HM- ( 231 ) 
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After a little algebra we turn Eq.(226) into a closed equation for S(j) 



"*(^)-(^) + (^,)- « 

Introducing the projection operators orthogonally to cj> and <f>* : 

Q = 1 — l^oX^ol (233) 
Q* = 1-1^X^51 (234) 

we have S± = QS and 

( H GP + gN Q\M 2 Q gN Q4>lQ* \ 

\ -gN Q Q*(j>fQ -[H GP + gN Q\<f> \ 2 QY ) ( > 

where the Gross-Pitaevskii Hamiltonian H GP is defined in Eq.(228). In general the 
operator C is diagonalizable. 

6.1.3 Spectral properties of £ and dynamical stability 

Consider an eigenvector of C with eigenvalue e k : 

<::)-(::) « 

The free evolution of this mode, that is for SU = , is given by the phase factor 
exp(—iekt/h) . This factor remains bounded in time provided that the imaginary part 
of €k is negative, which leads to the dynamical stability condition 

lm(e k ) < for all k. (237) 

One has the following three interesting spectral properties. 

1. is also an eigenvalue of C . 

2. ^ ^ j is also an eigenvector of C , with eigenvalue —e% . 

3. | Uk ) is an eigenvector of £) with eigenvalue e k . 
V ~Vk J 
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The last two of these three properties can be checked by direct substitution. They can 
be viewed more elegantly as a consequence of the symmetry properties: 

(238) 
(239) 

As we shall see this last property involving £) is useful to write C in diagonal form. 

The first of these properties is easy to prove when ^ is real. In this case the operator 
C is a real operator (as Q , U Q , , A are real) so that complex eigenvalues come by 
pairs of complex conjugates. When </5> is complex one can use the following mathematical 
fact: if e& is an eigenvalue of an arbitrary operator C , e* k is an eigenvalue of the 
operator & . Here C) differs from £ by a unitary transform (239) so that it has 
the same spectrum as C . This property of the spectrum is actually known in classical 
mechanics for a linearized Hamiltonian system, and the Gross-Pitaevskii equation can be 
viewed as a classical Hamiltonian equation for a continuous set of conjugate coordinates 
q = Re((/)),p = lm(0) . 

As a consequence of this first spectral property of C the dynamical stability condition 
Eq.(237) can be reformulated as 

Im(e fc ) = for all k (240) 

that is all the eigenvalues of C have to be real to have a dynamically stable condensate 
wavefunction. We assume that this property is satisfied in the remaining part of this 
subsection 6.1. 



6.1.4 Diagonalization of C 

As C is not a Hermitian operator the eigenbasis of C is not orthogonal (see the minus 
sign in the second line of C , due to Bose statistics; one does not have this sign in the 
BCS theory for interacting fermions). 

To write £ in diagonal form the knowledge of the eigenvectors is then a priori not 
sufficient. One generally proceeds as follows: 
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Reminder: Let M be a diagonalizable but not necessarily Hermitian operator. Then 
the diagonal form of M can be written as 

M = Y t m k \^)(rt\ (241) 

k 

where is a right eigenvector of M with eigenvalue m& : 

M\tf) = m k \itf) (242) 
and | is a left eigenvector of M with eigenvalue : 

{rt\M = mM\ (243) 

or equivalently 

MM) = ™*M)- (244) 
The normalization of the left and right eigenvectors is such that 

<VM> = S k , k ,. (245) 

is then called the adjoint vector of . 
We apply this reminder to C . We have already defined the right eigenvector: 

|Vf> = ( ) , eigenvalue of £ : e k . (246) 

From the third of the above spectral properties of C we can easily obtain the correspond- 
ing left eigenvector up to a normalization factor: 

l^ L ) = A4 ( ^ ) , eigenvalue of O : e* = e k . (247) 

The normalization condition Eq.(245) imposes 

= 1 = Af k [(^\u k ) - (v k \v k )]. (248) 

It is then natural to normalize the right eigenvectors in such a way that the quantity 
between square brackets is ±1 , leading to A4 = ±1 . 

We therefore group the eigenvectors of £ in three families: 
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• the + family, such that (u k \u k ) — (v k \v k ) — +1 

• the — family, such that (uk\uk) — (vk\vk) = —1 

• the family, such that (u k \u k ) — (v k \v k ) = . 

From the spectral property number 2 we see that there is a duality between the + 
family and the — family. Conventionally we will refer to eigenvectors of the + family 
as (u k ,v k ) (eigenvalue e k ) and the eigenvectors of the - family will be expressed as 
(vl,u%) (eigenvalue -e k ). 

Generally there are only the following two members in the family: 

c°) - (;)• « 

One can check that these two vectors are eigenvectors of C with the eigenvalue zero. [In 
the case of Cqp one finds in general only one zero-energy eigenmode, the missing one 
leads to the non-diagonalizability] . In general these two vectors span the whole family. 
As they are also eigenmodes of the operator & with the eigenvalue zero they are actually 
their own conjugate vectors! From Eq.(245) we then get the important property: 

(0o|wjfe) = (0o|^) = for all k in + family. (250) 

It is important to note that the + in the denomination " + family" refers a priori to 
the sign of (u k \u k ) — (v k \v k ) and not to the sign of e k ! 

6.1.5 General solution of the linearized problem 

We expand the unknown column vector of Eq.(232) onto eigenmodes of C . We assume 
that the only modes of the family are the ones of Eq.(249); these modes do not 
contribute to the expansion as (<^o|<^±) = ($5 1 ^01) = ■ We tnen S et tne expansion 



with the coefficients 

h(t) = ((u k \, -(v k \) ( ) = / d¥ [ul(r )54> ± (f,t) - vlirWAm . (252) 
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As the second component of the expanded column vector is the complex conjugate of the 
first component the amplitudes on the — family modes are the complex conjugates of 
the amplitudes bj~ on the + family modes, that is b k . 

Similarly one expands the source term of Eq.(232) on the eigenmodes of C . The 
components on the + family modes are given by 

**(*) = / d 3 r K(f)S ± (f, t) + vl{r)Sl{r, t)] . (253) 

Note the absence of - sign here, due to the fact that the second component of the source 
column vector is the opposite of the complex conjugate of the first component. Finally 
the projection of Eq.(232) on the eigenmodes of the + family gives the set of equations: 

ih^h(t) = € k b k (t) + s k (t) (254) 
which can be integrated including the initial condition 8<j>±_{to) = : 

h(t)= f~ t0 %s k {t-T)e- i ^l\ (255) 
Jo in 

6.1.6 Link between eigenmodes of C GP and eigenmodes of C 

The linear operators C and Cqp describe the same physical problem, so that one expects 
in particular that their spectra, which correspond to the linear response frequencies of the 
condensate, are the same. 

This expectation is confirmed by the simple result, that one can check by direct 
substitution: if (\u k p ),\v k p )) is an eigenvector of Cqp with the eigenvalue then 
(Q\u k p ), Q*\v k p )) is an eigenvector of C with the same eigenvalue . 

6.2 Examples of dynamical instabilities 

We consider simple stationary solutions of the Gross-Pitaevskii equations that are dynam- 
ically unstable, that is with non-real eigenfrequencies e k /% . The situations considered 
correspond to a gas trapped in a cubic box with periodic boundary conditions; analyti- 
cal calculations can then be performed. The conclusions remain qualitatively correct for 
harmonic traps. 
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6.2.1 Condensate in a box 



The atoms are trapped in a cubic box of size L , and we assume periodic boundary 
conditions. An obvious solution of the time independent Gross-Pitaevskii equation is 
then the plane wave with vanishing momentum, 



1 



Mr) = 



It has a chemical potential 



(256) 
(257) 



H = gN \(j)o\ 2 = p g 

where p = N /L 3 is the density of condensate atoms. 

To obtain the linear response frequencies of the condensates we calculate the spectrum 
of Cqp , this operator takes here the very simple form: 



C GP — 



I ^ A 



V 



-po9 



Po9 
2m 



/ 



(258) 



Using the translational invariance of this operator we seek its eigenvectors in the form of 
plane waves: 



vf P (r) 



L 3 / 2 



(259) 



Within the subspace of plane waves with wave vector k , C GP is represented by the 2x2 
non-Hermitian matrix: 



2L2 



C GP [k] = 



\ 



2m 
~Po9 



+ Po9 



Po9 
rtfk 2 



2m 



+ Pv9 



\ 



J 



(260) 



For k ^ this matrix can be diagonalized, giving one eigenvector of + family, with the 
eigenvalue 

W/W Nil/2 



2m \ 2m 



(261) 
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and one eigenvector of — family with the eigenvalue — 
family can be chosen as 



The eigenvector of the + 



2m 
h 2 k 2 n 

/ h 2 k 2 



1/4 



-1/4 



2m 



2/,2 



hfk 

V 2m 



+ 2po9 



(262) 



with the correct normalization U% — V? = 1 . 

k k 

The spectrum Eq.(261) is the so-called Bogoliubov spectrum, as it was first derived 
by Bogoliubov. Physically it is a very important result. In the limiting case of the ideal 
Bose gas (g = 0) the spectrum is the usual parabola; for g ^ the spectrum is very 
different, and this deserves a more detailed discussion 

• Bogoliubov spectrum for g > 

In the case of effective repulsive interactions the Bogoliubov spectrum strongly differs 
from the one of a free particle in the low momenta domain h 2 k 2 /2m <C 2p Q g as it scales 
linearly with k : 

(263) 



~ Tik 



Po9 
m 



This linear dispersion relation leads to a propagation of low energy excitations in the 
condensate in the form of sound waves with a sound velocity c s given by 



c s = 



duj^ _ 1 de^ _ 



dk h dk 

or equivalently by the relativistic type formula 



mc 2 = p g. 



Po9 
m 



(264) 



(265) 



Superfluidity is an important consequence of this linear behavior of the spectrum at 
low k , as shown by an argument due to Landau and that we explain briefly. Consider a 
particle (of mass M ) sent in the atomic gas with an initial velocity u . The motion of this 
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particle can be damped by interaction with the condensate only if the particle can create 

some excitation of the condensate. Imagine that such an excitation is produced, with 

— » —> 
momentum k ; the particle experiences a momentum recoil of —hk and conservation of 

energy imposes 



1 



er = =-Mu 2 - 



Mu-hk 



12 . -* hfk 



k 2 2M L 

The velocity u has then to satisfy the inequality: 



2/,2 



= hk • u — . (266) 

2M v } 



u > 



k • u 



k 



>J> C *- ( 267 ) 



So a particle with an incoming velocity smaller than the sound velocity can move through 
the condensate without damping, only scattering on thermal excitations of the gas can 
damp its motion. This prediction has received an experimental confirmation at MIT [41]. 

At high momenta ( h 2 k 2 /2m p Q g ) corresponding to a velocity hk/m much larger 
than the sound velocity the Bogoliubov spectrum reduces to a shifted parabola 

e ^^r +Po ^^r +2po ^- /i - (268) 

This approximate form can be obtained by a series expansion of the general formula 
Eq.(261). It can also be derived more instructively from the observation that the off- 
diagonal coupling p Q g between the component and the component in the 2x2 
matrix Eq.(260) becomes very non- resonant at high k (because the diagonal terms for 
and have opposite signs); neglecting this coupling one recovers Eq.(268) with 

This last high energy property applies also in a non-uniform trapping potential: ne- 
glecting the off-diagonal coupling between u k and v k one approximates the high energy 
part of the Bogoliubov spectrum by the eigenvalues of 

h 2 

-—A + U + 2^ |^ | 2 - p (269) 
2m 

which (up to the shift —jj, ) is the Hartree-Fock Hamiltonian for non-condensed particles 
Eq.(116) in a regime of an almost pure condensate where the density of non-condensed 
particles is negligible as compared to the condensate density A^ol^ol 2 ■ 

From this we expect that the Hartree Fock approach is invalid for the low energy 
fraction of the non-condensed gas (energy typically less than this may become a 
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problem at temperatures k B T < \i , where one has to use the more precise Bogoliubov 
approach of §7. 

• Case of a negative g 

In the case of effective attractive interactions between particles the dynamical stability 
condition Im(e^) = is satisfied if and only if 

H 2 k 2 

+ 2#p > for all A; > 0. (270) 

2m 

If one considers the thermodynamical limit of an infinite number of condensate atoms 
with a fixed mean density p Q = N /L 3 the stability condition cannot be satisfied as k 
can be arbitrarily close to in an infinite box. One may then be tempted to conclude 
that condensates with effective attractive interactions cannot be obtained experimentally, 
attractive interactions leading to a spatial collapse of the gas. 

Experiments with atomic gases can deal however with small number of atoms and 
the simplifying assumption of a thermodynamical limit is not necessarily a good approx- 
imation. In the cubix box of size L with periodic boundary conditions the compo- 
nents of the wavevector k of an atom are integer multiples of 27r/L . The smallest but 
non zero modulus of wavevector that can be achieved is therefore 27r/L (by taking e.g. 
k x = 27r/L, k y = k z = ). Dynamical stability condition Eq.(270) can then be rewritten 



or equivalently in terms of the scattering length as 



Condensates for a < can contain a limited number of atoms proportional to the size 
of the condensate. 

Condition Eq.(271) has a clear physical interpretation: the energy gap in the spectrum 
of a particle in the box between the ground state and the first excited states should 
be larger than the mean field energy per particle: stabilization against collapse is thus 
provided by the discrete spectrum of the atoms in the trapping potential. This condition 
can thus be qualitatively extended to the case of an isotropic harmonic trap, Hco > 



as 




(271) 
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\g\N /a^ o where u is the oscillation frequency of the atoms in the trap and a^o = 
(h/muj) 1 / 2 is the typical spatial extension of the ground state of the trap. One then 
recovers up to a numerical factor the results of §5.2.1. 

6.2.2 Demixing instability 

We consider here atoms with two internal states a, b ; this model is relevant for exper- 
iments performed at JILA on binary mixtures of 87 Rb condensates, and also (if one 
includes a third atomic internal level) experiments at MIT in Ketterle's group on spinor 
23 Na condensates. 

To describe the elastic interactions between the atoms with two internal states one 
needs three coupling constants, all positive in the case of 87 Rb: g aa and g bb for inter- 
actions between atoms in the same internal state, g ab for interactions between atoms in 
different internal states: 



gaa 

g^ 

gab 



a + a -¥ a + a 

b+b^b+b (273) 
a + b — > a + b. 



In the JILA experiment internal states a and b correspond to different hyperfine levels of 
the atoms so that inelastic collisions such as a+a — > a-\-b are either strongly endothermic 
(and do not take place) or strongly exothermic (and result in losses of atoms from the 
trap); we neglect these inelastic processes. 

Omitting for simplicity the regularizing operator in the pseudo-potential we write the 
interaction Hamiltonian between the atoms in second quantized form as 

tfint = / d*f [HfdKf )#(r )Ur )k{r ) + |^(f )$(f )^ 6 (f ) 

+ ^4 f (r)^(r)4(r)^(r)] (274) 

where $ a and ^ are the atomic field operators for atoms in state a and b respectively. 
Note the absence of factor 1/2 in the a — b interaction term, which is best understood 
in first quantization point of view: all the pairs of atoms of the form a : i,b : j , where 
% running from 1 to N a labels the atoms in a and j running from 1 to N b labels the 
atoms in b , are different so that there is no double counting of these interaction terms 
and no factor 1/2 is required. 
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Using the trick of §5.1.3 we can rapidly derive the Gross-Pitaevskii equations for the 
condensate wavefunctions (j) a in state a and (j) b in state b , both wavefunctions being 
normalized to unity. One simply has to write the Heisenberg equations of motion for the 
field operators and perform the substitution 



N l J 2 fa 
Nl /2 fa 



(275) 
(276) 



where N ajb are the number of particles in condensates a, b . As in §6.2.1 we restrict to 
the case of atoms trapped in a cubix box of size L with periodic boundary conditions. 
We obtain the coupled time dependent Gross-Pitaevskii equations 



ihd t (j) a = 
ihdtfa = 



-—A + N a9aa \(t> a \ 2 + N b9ah \(t> b \ 2 
' ft 2 

-^A + N b9bb \<t> b \ 2 + N agab \<t> a \ 2 



fa 
fa- 



(277) 



Consider now a steady state solution of these coupled Gross-Pitaevskii equations. As 
we have not introduced any coherent coupling between the internal states a and b (no 
i^l^a term in the Hamiltonian) (j) a and (j) b can have in steady state time dependent 
phase factors evolving with different frequencies: 



fa(r,t) = <l>*fl{?)e-W 
fa(r,t) = M (Oe- iMt/ft 



(278) 
(279) 



From a more thermodynamical perspective we can also observe that the number of parti- 
cles N a and N b are separately conserved by the three elastic interactions of Eq.(273) so 
that two distinct chemical potentials \i a and \i b are required to describe the equilibrium 
state. 

The time independent Gross-Pitaevskii equations for fa ;0 and fa$ in the box have 
the natural solutions 

0a,o(O = MO = ^ (280) 
leading to the following expressions for the chemical potentials: 



Pa — Pa,o9aa + Pb,o9bb 
pb — Pb,09bb + Pa,09ab 



(281) 
(282) 
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where p a;b are the condensate densities in a, b . Although we assume here that all the 
coupling constants are positive it is physically intuitive that these spatially uniform solu- 
tions should become instable when the interactions between a and b are very repulsive; 
one then feels that the two condensates a and b have a tendency to spatially separate. 

To test this expectation we linearize the time dependent coupled Gross-Pitaevskii 
equations around the steady state, setting 



We obtain 



Mr,t) = e-^ n [0 a;O (r) + ^ a (f^)] 
Mr,t) = e-^[^ (r) + ^(f,t)]. 



n 2 

ifld t S(j) a = -—A5(/> a + p a ,Q9aa[^(t>a + S<f>* a ] + Pb,09ab[$<f>b + Hb\ 



(283) 
(284) 



(285) 



and a similar equation for 8(j) b exchanging the role of a and b indices. We look for eigen- 

modes of these linear equations, with eigenfrequency e/% and a well defined wave vector 
— # 

k . This amounts to performing the substitutions 

and equivalent changes for the b components. This leads to the eigensystem 



2^21 



e — 



2m 



2^21 



—e — 



2m 



Va = Pa,09aa [u a + V a ] + pb,09ab [v>b + V b ] (286) 



and similar equations obtained by exchanging the indices a and b . Taking as new 
variables the sums and the differences between u and v and using the first identity in 
Eq.(286) we eliminate the differences as functions of the sums: 

2me 



n 2 k 2 



[u a + v a ] . 



(287) 



We get from the second equality in Eq.(286) (and a b ) a two by two system for the 
sums u a + v a , u b + v b : 

I 2m Wk 2 ) }\u b + v b 



= 



(288) 
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where we have introduced the two by two matrix 



M= Pa,09aa Pb,09ab j (289) 
\ Pa,09ab Pb,09bb J 



This leads to the following condition for the spectrum: 

W fh 2 k 2 



e 2 = 



2m \ 2m 



+ 2^,2 



(290) 



where 771,2 are the eigenvalues of M . 

In the thermodynamical limit the mixture of condensates with uniform densities is 
dynamically stable provided that both eigenvalues 771,2 are positive. This is equivalent 
to the requirement that the symmetric matrix M is positive. As g aa > here this is 
ensured provided that the determinant of M is positive: 

det M = p aj oPb,o [g a a9bb ~ gib] > 0- (291) 
The mixture of spatially uniform condensates is therefore stable if 

9ab < (9aa9bb) 1/2 - (292) 

In this case one can check that the spectrum of the + family is given by the positive 
solutions e toEq.(290). The spectrum of the binary mixtures of condensates is then made 
of two branches, which are both linear at low momenta with sound velocities (rji^/m) 1 / 2 . 

What happens when this stability condition is not satisfied ? The condensates a 
and b have a tendency to separate spatially. This happens in the JILA experiment [42]. 
Actually our model in a box is too crude to be applied to the experimental case of particles 
in a harmonic trap, the stability condition Eq.(292) being marginally satisfied for 87 Rb; 
a numerical solution of the time dependent Gross-Pitaevskii equations is required in this 
case [43]. 

The occurrence of demixing when Eq.(292) is violated can be connected with the 
following simple energy argument. Consider a demixed configuration with all the N a 
atoms in the left part of the box in a volume i/L 3 and all the N b atoms in the right part 
of the box in the complementary volume (1 — v)L 3 . The condensate densities vanish on 
a scale on the order of the healing length £ of the gas, this leads to "surface" kinetic and 
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interaction energies negligible in the thermodynamical limit as compared to the volume 
interaction energy 

N l9aa N b 2 g bb 

2vL* + 2(1-^)L3- [2 ^ } 
We minimize this energy over v to obtain 

^demix = ^3 [Nlgaa + g bb + 2N a N b (g aa9bb ) 1 / 2 ] . (294) 

We find that the demixed configuration has an energy lower than the one of the spatially 
uniform configuration 

£unif = ^ [Nlgaa + g bb + 2N a N b9ab ] (295) 
precisely when the stability condition Eq.(292) of the uniform configuration is violated. 



6.3 Linear response in the classical hydrodynamic approxima- 
tion 

We consider in this subsection the case of a condensate in a harmonic trap. The eigen- 
modes of the linearized Gross-Pitaevskii equation can then in general be determined only 
numerically. In the Thomas-Fermi regime however approximate results can be obtained 
for the low energy eigenmodes of the system from the classical hydrodynamic approach, 
as we shall see now. 



6.3.1 Linearized classical hydrodynamic equations 

The classical hydrodynamic equations for the position dependent condensate density p(f ) 
and velocity field v(f) have been derived in §5.3.4: 

d t p + div[pv] = (296) 
m (d t + v • grad) v = -grad[f7 + pg\. (297) 

We linearize these equations around their steady state solution with density p and 
vanishing velocity field v* = in the unperturbed trapping potential U Q . Writing the 
trap potential as a perturbation of Uq : 

U(r,t) = U Q (r) + SU(f,t) (298) 
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and splitting p and v as 



p(f,t) = p Q (r) + 5p(r,t) 
v(f,t) = Q + 5v(r,t) 



(299) 
(300) 



we obtain the linearized equations: 



d t Sp + div [poSv] = 
md t 5v + grad [5pg] = — grad 5U. 



(301) 



Taking the time derivative of the first equation we obtain a term d t 5v that we can 
eliminate with the second equation. This results in a closed equation for the perturbation 
of density: 



d 2 Sp - div 



Po9 
m 



grad dp 



= div 



^grad SU 
m 



(302) 



The homogeneous part of this equation can be rewritten in a more suggestive way by 
introducing the position dependent velocity c s given by 



mc 2 s (r) = po{r)g. 
The homogeneous part of Eq.(302) then reads 

d 2 Sp — div c 2 (f )grad Sp 



(303) 



(304) 



which corresponds to the propagation of sound waves with a position dependent sound 
velocity c s (f) . Note that the expression (303) could be expected from the result Eq.(265) 
obtained in the spatially homogeneous case. 

The propagation of sound waves in a cigar shaped trapped condensate has been ob- 
served in Ketterle's group at MIT; the condensate was excited mechanically by the dipole 
force induced by a far detuned laser beam focused at the center of the trap. It is instruc- 
tive to note that to predict theoretically the velocity of sound obtained in the experiment 
one has to carefully solve Eq.(302), rather than to take naively the sound velocity on the 
axis of the trap; the naive expectation would be wrong by a factor of \/2 [44, 45]. 



6.3.2 Validity condition of the linearized classical hydrodynamic equations 

The classical hydrodynamic equations have been obtained from the time dependent Gross- 
Pitaevskii equation in hydrodynamic point of view by neglecting the quantum pressure 
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term (see §5.3.4). If we keep this term and linearize the resulting equation we get extra 
terms with respect to Eq.(301). One of the extra terms is 

ti 2 

ASp (305) 



mpo 

which should be small as compared to the "classical" pressure term gSp : 

— \ASp\ < g\Sp\. (306) 

mpo 

If we denote by k a typical wavevector for the spatial variation of Sp , ASp ~ —k 2 Sp 
and the condition for neglecting the quantum pressure term reads 

^k 2 , ^ 

< 9Po- 307 

rn 

This condition can be rewritten in a variety of ways. It claims that the wavevector k 
should satisfy 

h£ < 1 (308) 

where f = h/(2mp Q g) 1 ^ 2 is the healing length of the condensate. Eq.(302) cannot be 
used to describe perturbations of the condensate at a length scale on the order of f or 
smaller. 

The validity condition can also be written as 

Hk <C mc s or fikc s <C mc 2 s = p Q g ~ p (309) 

In terms of the Bogoliubov spectrum for the homogeneous condensate this means that 
the wavevector k has to be in the linear part of the excitation spectrum. The energy of 
the corresponding eigenmode ~ hkc s has to be much smaller than p . Therefore only 
the eigenmodes of Eq.(302) with eigenenergy much less than p are relevant: 

e < p (310) 

the higher energy modes are not an acceptable approximation of the exact eigenmodes of 
the condensate. Note that as shown in [46, 47], Eq.(310) is a necessary but not sufficient 
condition in a trap. 
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6.3.3 Approximate spectrum in a harmonic trap 

We look for eigenmodes of the homogeneous part of Eq.(302) with eigenenergy e . They 
solve the eigenvalue equation 

-e 2 Sp = div [c 2 (r )gradSp] . (311) 

In the present Thomas-Fermi regime (? s is a quadratic function of the coordinates as 

it is proportional to the condensate density p Q . We can therefore solve the eigenvalue 

equation using an ansatz for Sp polynomial in the spatial coordinates x,y,z . If Sp is 

— » 

a polynomial of total degree n , grad Sp is a polynomial of total degree n — 1 ; after 
multiplication by (? s and action of the div operator we get a polynomial of total degree 
(n— l) + 2 — l = n. The subspace of polynomials of degree < n is therefore stable. 

For low values of the total degree n an analytical calculation is possible. For example 
in a general harmonic trap with atomic oscillation frequencies uj x ,uj y ,uj z one can check 
that the polynomials Sp = x,y,z (respectively) are eigenvectors with eigenvalues e = 
huj x ,fwj y ,huj z (respectively). These three modes correspond to the oscillation of the 
center of mass of the gas, which is exactly decoupled from all the relative coordinates of 
the particles in a harmonic trap; for this specific example the frequencies (but not the 
modes!) predicted by classical hydrodynamics are exact. These "sloshing" modes are 
used experimentally to determine accurately the trap frequencies w x ,y,z ■ 

Another important example is the case of a total degree n = 2 . One can check that 
the subclass of polynomials involving the monomials x 2 , y 2 , z 2 and 1 is stable, which 
corresponds to the ansatz 

Sp(f) = B+ £ A a rl (312) 

a=x,y,z 

By inserting this ansatz into Eq.(311) we arrive at the eigenvalue system 

(e/h) 2 A a = 2u 2 a A a + ulY,Ap. (313) 

This eigenvalue system can be obtained more directly as in Eq.(216) by a linearization 
of the equations Eq.(212) for the scaling parameters X a around their steady state value 
X a = 1 . Physically the modes identified by the ansatz (312) are therefore breathing 
modes of the condensate. The frequency of one of these modes (the one of §5.4.3) has 
been measured with high precision at MIT in a cigar-shaped trap, it differs from the 
Thomas- Fermi prediction ~ (5/2) 1 / 2 lu z (see Eq.(218)) by less than one percent [48]. 
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In the case of an isotropic harmonic trap all the eigenenergies e can be calculated 
analytically (keeping in mind that modes with e > \i are not properly described by 
classical hydrodynamics !). As shown in [49] one uses the rotational symmetry of the 
problem, as in the case of Schrodinger's equation for the hydrogen atom, with the ansatz 



where 9, (j> are the polar and azimuthal angles of spherical coordinates, Y™ is the spher- 
ical harmonic of angular momentum I . The last factor Pi, n {r) is a polynomial of degree 
n in r . As P^ n is a polynomial the recurrence relation obtained from Eq.(311) for 
the coefficients of the monomials r j should terminate at j = n . This leads to the 
eigenfrequencies f2 = e/% such that 



where uj is the oscillation frequency of the atoms in the trap. For I = l,n = we 
recover the sloshing modes with frequency = u . 



7 Bogoliubov approach and thermodynamical stabil- 
ity 



Imagine that we have already checked the dynamical stability of a steady state solution 
4>q of the Gross-Pit aevskii equation for the condensate wavefunction. We cannot relax 
yet and be sure that the condensate will remain in state </5> in the long run. 

What is missing is a check that interaction of the condensate with the non-condensed 
cloud does not induce an evolution of the condensate far from the predicted ^ ■ We 
have to check what is called thermodynamical stability of the condensate as it involves 
the "thermal", non-condensed component of the gas. 

This check will be performed in the low temperature domain ( T <C T c ) using the 
Bogoliubov approach. We will then present examples of thermodynamical instabilities. 

We give here a summarized account of the U(l) -symmetry preserving Bogoliubov ap- 
proach developed in [50, 51]. In contrast to the almost general attitude in the literature 
this approach does not assume a symmetry breaking state with (t/>) ^ but consid- 
ers instead a state with a fixed total number of particles. A different U(l) -symmetry 
preserving Bogoliubov approach was developed long ago in [52]. 



6p(f) = YT(0, <t>yPi, n (r) 



(314) 



n 2 = (2n 2 + 2nl + 3n + l)io 2 for any n, I > 



(315) 
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This allows to eliminate a technical problem of the symmetry breaking approach: if 
{4>)(t = 0) ^ the state of the system necessarily involves a coherent superposition of 
states with different total number of particles; such a state cannot be stationary (as states 
with different number of particles have also different energies) and it experiences a phase 
collapse — > making the description of the evolution of the system more involved. 

7.1 Small parameter of the theory 

We restrict in this section to a steady state regime where most of the atoms of the gas 
are in the condensate. We split the atomic field operator as 

^(f) = ^(f)fi^ + <ty(f) (316) 

where 4> Q is the condensate wavefunction and a<p annihilates a particle in the mode </5> . 
The idea is to treat Sijj as a perturbation with respect to 0o^ o : Let us compare indeed 
the typical matrix elements of these two operators: 

6$~ (6$ 6$) 1/2 ~ (p') 1/2 (317) 

where p' is the density of non-condensed particles whereas 

~ Wo 1/2 0o ~ pj /2 (318) 

where p is the condensate density. We will therefore assume 

p' < Po (319) 

and even more 

N' = J d*rp'{r) < N ~ N (320) 

where N is the total number of particles in the gas. Using these two assumptions a 
systematic expansion of the field equations in powers of the small parameter 

e = (N'/No) 1 ' 2 < 1 (321) 

can be performed. We give here a somewhat simplified presentation. 
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The following identities are important properties of Sip . First, 6$ is the part of the 
atomic field operator transverse to the condensate wavefunction: 

J d 3 r(/>o*(r)5ip(r) = 0. (322) 

As a consequence the bosonic commutation relation obeyed by Sijj involves matrix ele- 
ments of the projector Q orthogonal to 4> Q rather than the identity operator: 

), Sft(r 2 )] = (n \Q\f 2 ) = S(n - r 2 ) - M(? 2 ). (323) 

Second, there should be no coherence in the one-body density matrix between the con- 
densate and the non-condensed modes, or equivalently </5> should be an eigenstate of the 
one-body density matrix (with eigenvalue N ): 

= 0. (324) 

This last identity is used to calculate the condensate wavefunction order by order in the 
small parameter [50, 51]. 



7.2 Zeroth order in e : Gross-Pitaevskii equation 

To zeroth order in e all the particles are supposed to be in the condensate. In steady 
state one then recovers the time independent Gross-Pitaevskii equation with the further 
approximation N Q ~ N : 

h 2 



-^A + t/ + ^|0o| 2 



^o- (325) 



7.3 Next order in e : linear dynamics of non-condensed particles 

Calculation to first order in e corresponds to a linearization of the Heisenberg field equa- 
tions around (pod^ keeping terms up to first order in Sijj . Equivalently it corresponds 
to a quadrat izat ion of the Hamiltonian around 0o^ o keeping terms up to second order 
in 5ip . 

We use this quadratization approach here. At this order of the calculation the reg- 
ularizing operator of the pseudo-potential can be neglected, divergences due the use of 



Bogoliubov approach 



89 



the non-regularized gS potential coming at the next order e 2 . We therefore take the 
Hamiltonian 

H = Jd 3 f $hii> + |^t^, . (326) 

The one-body Hamiltonian hi contains the kinetic energy and the trapping potential 
energy: 

h 2 

hi = -—A + U. (327) 
2m 

It does not contain any —\i term as we use here in the canonical rather than grand 
canonical point of view, the total number of particles being fixed to N . 

We substitute expansion (316) for Eq.(326) and we keep terms up to quadratic in Sip . 

• The contribution of hi is quadratic in $ so that all the terms should be kept. One 
of the contributions is 

(| d 3 f^ ^i^o)4o^o- (328) 

One can use the following trick to express this quantity in terms of Sip : from Eq.(322) 
one sees that the total number of particles operator can be written as 

N = J d 3 rftip = h<t> + SN (329) 

that is the sum of the number operator of condensate particles: 

= 4o^o ( 33 °) 

and the number operator of non-condensed particles: 

SN = J d 3 rS$Sip. (331) 

We therefore obtain 

= N - SN. (332) 

• The expansion of the interaction term gives the following up to quadratic contributions: 

k 

•fit fit 

+ 2#V^o4o^ t( ^^o- ( 333 ) 



+ g[(p *(po*(podl o a^ o a^ Sip + h.c] 
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The first line of this expression is transformed using Eq.(332): 

4o4o Wo = - 1) = & (ft - 1) - 2 ft 6 ft + ■ ■ ■ ( 334 ) 

with N ~ N — 1 . When we group the above term in SN with the one coming from hi 
and we replace ./V by TV we obtain 

—SN [| d 3 f o *Mo + Ng\</>o | 4 = —fiSN (335) 

as 0o solves the Gross-Pitaevskii equation (325). This is amusing: we obtain formally a 
grand canonical Hamiltonian for the non-condensed particles, the reservoir being formed 
by the condensate particles! In the second line of Eq.(333) we replace a^a^ by N as 
the corrective term SN would lead to a cubic contribution in Sip . 

Another important transformation is performed by collecting the terms linear in Sip 
from Eq.(333) and from the contribution of hi , leading to 

f d 3 rct> *al o [hi + g\c/>o\ 2 N}S^ + h.c. (336) 

As (j) Q solves the Gross-Pitaevskii equation, the operator between brackets, when acting 
on the left on <^ * , gives a contribution /^ * orthogonal to Sip (see Eq.(322)). The 
sum of all the terms linear in Sip therefore vanishes! We shall see later that this could 
be expected from Eq.(324). 

7.4 Bogoliubov Hamiltonian 

We now collect all the terms of H up to quadratic in Sip and express them in terms of 
the field operator 

A ( f ) = (337) 

This operator commutes with the operator N giving the total number of particles: it 
transfers one non-condensed particle into the condensate. As the operator Sip , it is 
transverse to (p Q (see Eq.(322)). By definition of the condensate wavefunction A has a 
zero mean (see Eq.(324)). 

In general it is difficult to exactly eliminate Sip in terms of A . Fortunately at the 
present order of the calculation we can use the assumption of a very small non-condensed 
fraction so that one has for example 

6$ 6$ ~ S$a^ N- l a\ Q Si) = A+A . (338) 
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To the same order of approximation A obeys the same commutation relation Eq.(323) 
as S'lp . 

The final result can be written in term of the operator C introduced in §6, with the 
approximation N Q ~ N : 

tfquad = f(N) + lf cPr (At, -k)C ( £ j (339) 

where the function / is specified in [51]. 

From this quadratic Hamiltonian the Heisenberg equations of motion for the field A 
have the suggestive form 

A Ma.)- « 

We note that an hypothetic term of Hq Ua d linear in A would give rise to a source term 
in Eq.(340) preventing one from satisfying Eq.(324) at all times! 

The result Eq.(340) is really a crucial one. It shows that the linearized evolution of the 
non-condensed part Sijj of the atomic field is formally equivalent to the linearized response 
of the condensate to a classical perturbation (e.g. of the trapping potential) derived from 
the Gross- Pit aevskii equation; both treatments indeed exhibit the same operator C (see 
Eq.(232)). 

All the machinery of §6 can therefore be used. We expand the field operator A on 
the eigenmodes of C . We assume here dynamical stability and that the only eigenmodes 
in the family are the zero-energy modes (<^o,0) and (O,0 O *) to which the field A is 
orthogonal according to Eq.(322). We therefore get an expansion similar to Eq.(251): 

with the important difference that the coefficients in the expansion are now operators: 

b k = J d 3 f [ul(f)A(r ) - vl(f)A\f)\ . (342) 

From the normalization condition between the eigenvectors of £ and their adjoint vectors 
(see §6) one shows that the b k obey bosonic commutation relations: 

ML] = *M' and M*'] = 0. (343) 
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bk corresponds formally to an annihilation operator; as \vk) ^ in general b k does 
not simply annihilate a particle as it is a coherent superposition of A (which transfers 
one non-condensed particle to the condensate) and of A* (which transfers one condensate 
particle to the non-condensed fraction) . One then says that b k annihilates a quasi-particle 
in mode k . 

Finally we rewrite the Hamiltonian Eq.(339) in terms of the b k 's: 

H vuA = E {N) + £ e k b\b k . (344) 

fcG+family 

We recall that e& is the eigenenergy of the mode (uk,Vk) of the + family. Our quadratic 
Hamiltonian describes a gas of non-interacting quasi-particles: this is the so-called Bo- 
goliubov Hamiltonian. 

The ground state of H qua d is obtained when no quasi-particle is present, it corresponds 
to all the modes k being in vacuum state: 

£ q uad|0) = E o (N)\0) (345) 

where 

Sjfc|0> = V/c. (346) 

E is therefore the Bogoliubov approximation for the ground state energy of the gas. 
To get a finite expression for E one has to include the regularizing operator in the 
pseudo-potential . 

The excited states of the system are obtained in the Bogoliubov approximation by 
successive actions of the b\ 's. For this reason b\ is said to create an elem en tary excitation 
k in the system, to distinguish with collective excitations involving all the particles of the 
condensate (induced e.g. by a perturbation SU of the trapping potential). We emphasize 
again that the elementary excitations of the gas have the same frequency e k /h as the 
collective excitations in the linear response domain ( SU small enough). This intriguing 
property is valid only at the presently considered regime of an almost pure condensate 
(T<T C ). 

7.5 Order e 2 : corrections to the Gross-Pitaevskii equation 

Expanding the Heisenberg field equations for ip keeping terms up to N l f 2 e 2 one can 
calculate the first correction to the prediction Eq.(325) for the condensate wavefunction. 
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This correction includes (i) the fact that the number of condensate particles N Q rather 
than the total number of particles N should appear in the Gross-Pitaevskii equation, 
and (ii) the mechanical back-action of the non-condensed particles on the condensate in 
the form of mean field potentials. 

The calculations are a bit involved [51] and require the use of the regularizing operator 
in the pseudo-potential (a fact realized in [53] but not yet in [51]). We give here only the 
result. The condensate wavefunction is given by an expansion 

O = O (O) + O + (e 2 ) (347) 

where , zeroth order approximation in e , is the solution of Eq.(325). The correction 
0o ^ is of order e 2 ; its component on is purely imaginary (as both O and 

are normalized to unity) and can be considered as a (not physically relevant) change of 
global phase of ; the part of orthogonal to is given by: 



(2) 



- £ u*^rl-3^J (348) 

where Q projects orthogonally to O and where the source term S is equal to 

5(f) = " (l + (SN))g\</> W(r)W(n 

+ 2g{k\r )A(r ))<j>^\?) + 9 [d* (*(A(r - s/2)W+ s/2)))] s ^<t>o {0> (f) 

- g Jd 3 s\4> {0) (?)\ 2 ([tt(s)cf> W(s) + A(s)V 0) *(s)] A(f )}. (349) 

The first line of this expression contains the effect of the depletion of the condensate, 
the number of non-condensed particles (SN) being calculated in the Bogoliubov approx- 
imation, see discussion in §7.7. The other terms are mean field terms, among which one 
recognizes the Hartree-Fock contribution 2g{k\r )A(f )) already obtained in §4. 



7.6 Thermal equilibrium of the gas of quasi-particles 

In the Bogoliubov approximation the quasi-particles behave as an ideal Bose gas; such a 
gas can reach thermal equilibrium only by contact with a thermostat. There is no such 
thermostat in the experiments on trapped Bose gas, relaxation to thermal equilibrium 
has to be provided instead by interactions between the atoms. 
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Fortunately the full Hamiltonian Eq.(326) contains terms cubic and quartic in Sip : 
when expressed in terms of the b 's and ft 's they correspond to interactions between the 
quasi-particles which will provide thermalization. Two situations can then be considered, 
depending on the sign of e k . 

• "Good" case: e k > for all k in + family. We assume for simplicity thermal 
equilibrium in the canonical point of view (which should be equivalent to the micro- 
canonical point of view in the limit of large number of particles) with a N -body density 
matrix 

Pi n = ^e-<> A e^Vad . (350) 

^ -^quad 

We suppose therefore that the interactions between quasi-particles, essential to ensure 
thermalization, have a weak effect on the thermal equilibrium state. From Eq.(350) we 
finally obtain the mean number of quasi-particles in mode k : 

= iSTTT ( 351 ) 
This good case corresponds to a thermodynamically stable condensate in state 4> Q . 

• "Bad" case: there is a mode k Q in the + family such that e ko < . In this case the 
quadratic Hamiltonian ff qua d contains a harmonic oscillator of frequency uj ko = \ e k \/h 
having formally a negative mass M : 

-ho\blk = ^[P 2 ko +<Ql] (352) 

where Pk and Qk correspond formally to a momentum and position operator. By 
collisions with the quasi-particles of positive energy the mode /c can loose energy which 
increases its own excitation; if the number of quanta in the mode can become comparable 
to N the process of thermalization of the gas may lead the condensate to a state different 
from the predicted (f>o . 

This phenomenon of thermodynamical instability should not be confused with dynam- 
ical instability of the condensate, where e ko is complex; e.g. the case of a purely imaginary 
e ko corresponds formally to an oscillator in an expelling potential, [P ko — \oj ko \ 2 Ql Q ]/(2M) . 

7.7 Condensate depletion and the small parameter (pa 3 ) 1 / 2 

We assume the situation of thermodynamical stability. We calculate the mean number 
(SN) of particles out of the condensate, that is in all the modes orthogonal to the con- 



Bogoliubov approach 



95 



densate wavefunction (f> : 

(SN) = J d 3 r(5ft(r)5iP{r)). (353) 

In this way we can calculate explicitly the small parameter of the present theory given in 
Eq.(321). 

To lowest order in the Bogoliubov approximation we can replace Sip by A in the 
above expression: 

(SN) ~ J d 3 f (At(r )A(r)>. (354) 

We replace A by its modal expansion; using the approximation in Eq.(350) the only 
terms with non-zero mean are (b\b k ) given by Eq.(351) and 

(hbl) = {b\b k ) + 1. (355) 

We therefore obtain: 

{SN) * £ (blh){(u k \u k ) + (v h \v k )] 

A;G+family 

+ £ («*l«*>- ( 356 ) 

fcG+family 

The contribution of the occupation numbers (b\b k ) corresponds to the so-called thermal 
depletion of the condensate, as it is non-zero only at finite temperature. The contribution 
of the +1 coming from the identity (355) is the so-called quantum depletion of the 
condensate: it expresses the fact that, even at zero temperature, there is a finite number 
of particles out of the condensate due to atomic interactions, contrarily to the ideal Bose 
gas case where all the v k 's vanish. 

It is instructive to calculate the depletion of the condensate in the homogeneous case 
in the thermodynamical limit, replacing the sum over the wavevector k characterizing 
each mode of the + family by an integral. From Eq.(262) we calculate 

^ + pq 

(u k \u k ) + (v k \v k ) = 1 + 2(v k \v k ) = 2 ™ ™ (357) 

\t [\t + 2 HJ 

At zero temperature we obtained for the non-condensed fraction of particles the following 
integral: 



(SN) , N 16 / ox 1/2 r+™ i 2 T <? 2 + l/2 . . 

; ( r = °) = ^(^ 3 ) I 99 [^fv^ _ 1 (358) 



N 
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where we have integrated over the solid angle in spherical coordinates and we have made 
the change of variable 

This integral can be calculated exactly: 

fV = «) = ^ (360) 

In this way we obtain a very important small parameter of the theory, (pa 3 ) 1 ^ 2 , which 
characterizes the domain of a weakly interacting Bose gas. This small parameter is similar 
to the condition obtained in Eq.(95) of §3 already with a totally different point of view, 
a condition of Born approximation on the pseudo-potential! In the typical experimental 
conditions of MIT for sodium atoms we have p ~ 10 14 cm ~ 3 and a ~ 25 A, which leads 
to a small parameter (pa 3 ) 1 ^ 2 ~ 10 -3 . The result (360) also gives us the opportunity to 
recall that the number of non-condensed particles should not be confused with the number 
of quasi-particles for an interacting Bose gas: one notes here that the second quantity 
vanishes at T = whereas the first one does not. 

At finite temperature there is, in addition to the quantum depletion, a thermal deple- 
tion of the condensate: 

<f> m - <f><«, . £f « &±p [ex P (jW + „*) - ,] - . ,3a:, 

The integral over q depends only on the parameter pg/iksT) . At low temperature 
UbT <C pg this parameter is large so that the modes with a q ~ 1 have a very low 
occupation number: one can neglect q as compared to one and one-half (which amounts 
to restricting to the linear part of the Bogoliubov spectrum) and one obtains a small 
correction to the zero temperature case: 



<^> (T) = M> (o) 




(362) 



At high temperature k B T ^> pg the major fraction of the populated Bogoliubov modes 
have a q much larger than unity, so that we can now neglect one and one-half as compared 
to q 2 (which amounts to restricting to the quadratic part of the Bogoliubov spectrum). 
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To this approximation we recover the ideal Bose gas result 



dB 



where AdB is the thermal de Broglie wavelength and £ stands for the Riemann Zeta 
function (see §2). For a given temperature indeed the density of non-condensed particles 
in presence of a condensate saturates to its maximal value C(3/2)A^"| for the ideal Bose 
gas in a box (see Eq.(5) and Eq.(2)). 

In the high temperature regime ksT ^> pg our Bogoliubov approach can therefore be 
valid only if pA^ B > C(3/2) . Both inequalities on p can be satisfied simultaneously only 
if 2((3/2)a/A dB <C 1 . Amusingly this condition is similar to the condition a A A; <C 1 
obtained in §3 (see Eq.(96)) for the validity condition of the Born approximation for the 
pseudo-potential . 



7.8 Fluctuations in the number of condensate particles 



The Bogoliubov theory that we have presented also allows a calculation of the fluctuations 
of No in the canonical ensemble. This has the advantage of removing the effect of 
fluctuations in the total number of particles present in the grand-canonical ensemble, a 
"trivial" contribution to the fluctuations of N Q . 

As the total number of particles is fixed to N the variance of N is exactly equal 
to the variance of SN , number of non-condensed particles. The mean value of SN has 
been given in the previous section, we now have to calculate the mean value of its square: 

((SN) 2 ) = J 'dnj 'dfa^^)^^!)^^)^^)) (364) 
= (SN) + J dn J df 2 (6ft {n )5$(r 2 )<fy(fi )<ty(f 2 )) (365) 

where we have used the commutation relation Eq.(323) and the fact Eq.(322) that (j) is 
orthogonal to Sip . To lowest order in the Bogoliubov approximation we can replace Sip 
by A in the above expression and take the approximation Eq.(350) for the equilibrium 
density operator. As -£/q U ad is quadratic in the field A , Wick's theorem can be used. We 
finally find for the variance of the number of non-condensed particles: 



Var(cW) ~ (SN) + J ' dfi j df 2 \(tt(n )A(f 2 )>| 2 + |(A(fl )A(f 2 )) 



(366) 
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The variance in Eq.(366) is simple to calculate in the homogeneous case of a gas in 
a cubic box of size L with periodic boundary conditions. The eigenmodes (u, v) are 
simply plane waves Eq.(259) with real coefficients U k , V k given in Eq.(262) and depending 
only on the modulus k of the wavevector k . We find for the two correlation functions 
appearing in Eq.(366) the simple expression: 

<At(fl)A(? 2 )> = ^E[(^ + ^H + ^]e* (ft - ft) (367) 
(A(rl )A(f 2 )) = ij;£4V,(l + 2n,)^- s ) (368) 

where we have introduced the mean occupation numbers n k = (blb k ) . After spatial 
integration of the square modulus of these quantities we obtain 

Var(<W) = (5N) + £ [{Ut + if) n k + vff + U 2 k V*(l + 2n k f (369) 

where the mean number of non-condensed particles (SN) is already given in Eq.(356): 

(5N) = Y,{UZ + V k 2 )n k + V*. (370) 

We first apply formula (369) to the limiting case of zero temperature. All the occu- 
pation numbers n k vanish. For the ideal Bose gas ( g = ) all the V k 's are zero and 
the variance of SN vanishes as expected, since all the particles are in the ground state 
of the box. For the interacting Bose gas we find 

Va r (^)(T = 0)a|:^-L_ (371 ) 

where q is given as function of k by Eq. (359) . In the thermodynamical limit we replace 
the discrete sum by an integral and this leads to a variance scaling as the number of 
particles for a fixed density: 

Va^A>)(T = 0) =27rl/2(pa3)1 ^ (372) 

In the regime of validity of the Bogoliubov approach one has pa 3 < 1 so that the 
fluctuations of AT are sub-poissonian. 
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The situation can be totally different at finite temperature. Consider first the case of 
the ideal Bose gas [54]: 

V*r(6N) = 5>*(l + »*) 

= g4sinhie,/2) (3?3) 

with 6k — h 2 k 2 / (2m) . In the thermodynamical limit one may be tempted to replace 
in the usual manner the sum over k by an integral. This leads however to an integral 
divergent in k — : the integrand scales as ljk A , which is not compensated by the 
Jacobian k 2 of three-dimensional integration in spherical coordinates. In this case the 
contribution of the sum in the thermodynamical limit is dominated by the terms close to 
k = where /5e^ <C 1 and the function sinh can be linearized. We then obtain 



Var(<W) 9=0 ^ V^. (374) 




In this expression we have introduced the kinetic energy difference between the ground 
state and the first excited state for a single particle in the box: 

and the sum ranges over all the vectors ft with integer components and a non- vanishing 
norm n . By a numerical calculation we obtain 

E^4 = 16-53 . . . (376) 

-t lb 

A remarkable feature is that the variance of N scales as L 4 that is as the volume of the 
box to the power 4/3 , or equivalently as the number of particles to the power 4/3 in 
the thermodynamical limit. This is much larger than N in the thermodynamical limit. 
Do the fluctuations of remain large in presence of interactions ? As the spectrum 
is linear at small k the divergence of n\ is only as 1/k 2 , which has a finite integral in 
three dimensions. However the mode functions £4, V& are also diverging in k — , each 
as 1/A; 1 / 2 , so that one recovers the 1/k 4 dependence of the summand close to k = . 
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As in the ideal Bose gas case we replace in the thermodynamical limit the summand by 
its low k approximation: 

* ~ m (377) 
u " ~ (378) 

and we keep in the summation the most diverging terms. We finally obtain [55, 56] 

Var(5jV) 9>0 4(^) 2 £l (380) 

Remarkably this result differs from the ideal Bose gas case Eq.(374) by a factor 1/2 only: 
fluctuations of N remain large. The fact that Eq.(380) does not depend on the strength 
g of the interaction is valid only at the thermodynamical limit and indicates that the 
limit g — > and the thermodynamical limit do not commute. 

7.9 A simple reformulation of thermodynamical stability condi- 
tion 

The thermodynamical stability condition simply requires that the Bogoliubov Hamilto- 
nian ff qua d is the sum of a constant (function of N ) and of a positive quadratic operator 
in the field variables. In the diagonal form (344) positivity is clearly equivalent to the 
requirement e k positive for all k in the + family. How to express this condition from 
the non-diagonal form (339)? We simply have to rewrite Eq.(339) as 

ff q uad = f(N) + \J d 3 r-(At, A) V C ( £ ) (381) 

where rj is the operator 

,-(; .;) m 

so that rjC is an Hermitian operator. 
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The thermodynamical stability condition is therefore equivalent to the requirement 
that rjC is positive: 

tjjC > 0. (383) 

More precisely, as A is orthogonal to O , f]C has to be strictly positive in the subspace 
orthogonal to (|0 O ),O) and (0, |^ *)) ■ 

We can give a simple physical interpretation of this condition: (j) Q has to be a local 
minimum of the Gross-Pitaevskii energy functional 

r h 2 



E[4>, 4>*] = N I d 3 f 



2 Jgrad^| 2 + U{f)W)\ 2 + ^NgWW 



(384) 



which is the expression of §5 with the approximation N ~ N . Let us consider indeed 
the variation 8E of E from E = E[4> ,(j)Q*] up to second order in a small deviation 
Scj) of (j) from $q . 

The terms linear in 8(j) are given by: 



SEW = N J 



d 3 f 



2m 



grad 0o* • grad 8<f> + U (r )^o*^0 + Ngfa* fioScj) + c.c. 



(385) 



By integration by parts and using the fact that </5> solves the Gross-Pitaevskii equation 
Eq.(325) we rewrite this expression as 

SEW = Nii[((t) Q \5(t)) + (S^o)]. (386) 

As both cf> and 0o are normalized to unity 8(j) actually fulfills the identity 

M,|ty> + <ty|to) = -(<W$. ( 387 ) 

The a priori first order energy change is a posteriori of second order: 

6E& = -Nii{6<l>\6<i>)\ (388) 

The terms a priori quadratic in 8<p are given by: 

r *2 



6E® = N f 



d 3 r 



2m 



grad 8<f)* • grad 8(j) + U{r)8<f)*8<f) + 2Ng\fa\ 2 8<i>*8<t> 



+ l -Ng^ 2 8cj> 2 + l -Ngcj) Q 2 8^ 2 



(389) 
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We transform this expression by splitting S(j) in a part parallel to O and a part orthog- 
onal to 4>q : 

5<j)(r) = 7 ^o(r ) + 8^{r ). (390) 

Using integration by parts, the Gross-Pitaevskii equation and the fact that operator C 
contains projectors orthogonally to <^o,0o* we are able to write SE^ as 

5EW = \^ Nfl + i( 7 + 7*) 2 jV 2 5 J d 3 r |^o| 4 + (7 + 7*)^ 2 <? / ^ItfoftV^x + ex.) 

+^iV / d¥ (ift, + M Id) ^ ) . (391) 

Note that we had to add /i times the identity matrix Id to tjjC as Eq.(389), contrarily 
to rjC , does not contain any term proportional to /i . From Eq.(387) we see that 7 + 7* 
is actually of second order in S(/)± so that it can be set to zero in Eq.(391). 

Summing the a priori first and second order energy changes we see that 8E^ exactly 
cancels the terms involving explicitly \i in Eq.(391) so that we arrive at 

SE =* ±nJ d 3 f(6<j>* ± , StJvC ( ) ■ (392) 

Thermodynamical stability that is positivity of tjjC is therefore equivalent to the Gross- 
Pitaevskii energy functional having a local minimum in </5> . 

7.10 Thermodynamical stability implies dynamical stability 

As we show now the positivity of rjC automatically leads to a purely real spectrum for C , 
that is to dynamical stability. Consider an eigenvector (u, v) of C with the eigenvalue 
e. Contracting the operator tjjC between the ket (\u),\v)) and the bra ((u\,(v\) we get 

({u\,(v\)vc( l ^=e[{u\u)-{v\v}]. (393) 

The matrix element of rjC is real positive as rjC is supposed to be a positive hermitian 
operator. We now face two possible cases for the real quantity (u\u) — (v\v) : 

• (u\u) — (v\v) = . In this case rjC has a vanishing expectation value in (\u), \v)) ; 
as rjC is positive, (\u), \v)) has to be an eigenvector of rjC with the eigenvalue 
zero; as rj is invert ible we find that (\u),\v)) is an eigenvalue of £ with the 
eigenvalue 0, so that e = is a real number. 
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(u\u) — (v\v) > : we get e as the ratio of two real numbers, so that e is real. 



7.11 Examples of thermodynamical instability 
7.11.1 Real condensate wavefunction with a node 

Let us assume that the solution of the Gross-Pitaevskii equation is a real function (/>o(r) . 
To decide if this solution is thermodynamically stable one has to check the positivity of 
the operator rjC . Consider an eigenvector of rjC with eigenvalue e : 



u 



= £ 




(394) 



This e should not be confused with the quasi-particle energies as r)C and £ have 
different spectra. By performing the sum and the difference of the two lines of Eq.(394) 
we get decoupled equations for the sum ip s = u + v and the difference ip d — u — v : 



2m 
2m 



+ U(r) + Ng<g{r) + 2NgQ^Hf)Q - }i 



+ U{?) + Ngtj>l{?)-n 



(395) 
(396) 



Both operators involved in these equations have to be positive to achieve thermodynamical 
stability. Note that for g > the positivity of the second operator Eq.(396) implies the 
positivity of the first one Eq.(395) as gQ(j>l{r)Q is positive. 

We therefore concentrate on Eq.(396). It involves the Gross-Pitaevskii Hamiltonian 



U GP = 1 ^ + U{r) + Ng4>l{r)- l i. 



(397) 



An obvious eigenvector of this Hamiltonian is ipd = 0o with eigenvalue e = , as 0o 
solves the Gross-Pitaevskii equation! The condition of a positive e in Eq.(396) simply 
means that O should be the ground state of Hop ! 

We can then invoke a theorem claiming that the ground state of a potential has no 
node [58]. If 4>o(r) has a node it cannot be the ground state of Hqp ■ The ground state 
of %gp has therefore an eigenenergy e lower than the one of O , that is lower than 
zero, so that the operator rjC is not positive and there is no thermodynamical stability. 
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As an example consider in a harmonic trap with eigenaxis z , a solution of the Gross- 
Pit aevskii equation even along x and y but odd along z , so that it vanishes in the plane 
z — . Such a solution exists, for g > : within the class of real functions (j) odd along 
z and even along x , y , the Gross- Pit aevskii energy E[<f>, (/>] , bounded from below, has 
a minimum, reached in cj) , and this O then solves the Gross-Pit aevskii equation. This 
solution however is no longer a local minimum of E[$,<f>*] when one includes all possible 
deviations of from <^ (complex and with no well defined parity along z ). 

7.11.2 Condensate with a vortex 

Can we get a thermodynamically stable condensate wavefunction with a node? To beat 
the results of the previous subsection we now assume <j> to be complex. 

A particular class of complex wavefunctions with a node are condensate wavefunctions 
with vortices. A vortex is characterized (i) by a nodal, not necessarily straight, line in 
(j) Q (the center of the so-called vortex core) and (ii) by the fact that the phase of O 
changes by 2qn , q non-zero integer, along a closed path around the vortex core ( q is 
the so-called charge of the vortex). This second property means that the circulation of 
the local velocity field (defined in §5.3.3) around the vortex core is 27rhq/m . 

A condensate wavefunction can have several vortices; the change of the phase of O 
along a closed path is now 27tq sum where q sum is the algebraic sum of the charges of the 
vortex lines enclosed by the path. 

It has been shown that a condensate wavefunction with a vortex in a harmonic trap 
is not thermodynamically stable [59]. In the limit of vanishing interaction between the 
particles (g = 0) this is clear indeed. Suppose that the trap is cylindrically symmet- 
ric with respect to z . |0o) can be chosen as (\n x = l,n y = 0,n z = 0) + i\n x = 
0,n y = l,n z = 0))/\/2 where \n x , n y , n z ) is the eigenstate of the harmonic oscillator 
with quantum number n a along axis a (a = x,y,z). The chemical potential is simply 
2huj X} y + ^huj z where uj a is the atomic oscillation frequency along axis a . One then 
finds that (\u) = \n x = 0,n y = 0,n z = 0), \v) = 0) is an eigenvector of rjC with the 
strictly negative energy e — —%uj Xiy . 

What happens in the opposite Thomas-Fermi regime of strong interactions? An intu- 
itive answer can be obtained in a 2D model of the Gross-Pitaevskii equation, assuming 
for simplicity a quasi-isotropic trapping potential and restricting to the following class of 
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condensate wavefunctions: 

<l>o(z,y) = &iow(z,2/) tanh[«|F - aR\]e ie ™. (398) 

In this ansatz (f) s \ ow (x, y) is the usual square root of inverted parabola Thomas- Fermi ap- 
proximation for a condensate wavefunction without vortex, with a radius R ; the tanh[ ] 
represents the correction to the modulus of 4>o due to the vortex core (of adjustable 
position aR and inverse width k ); 6$ R is the polar angle of a system of Cartesian co- 
ordinates (X, Y) centered on the vortex core, and represents (approximately for a ^ ) 
the phase of the unit-charge vortex. 

One then calculates the mean energy of ^> , with the simplification that 4> s \ ow (x,y) 
varies very slow at the scale of k~ 1 , and one minimizes this energy over n . This leads to 
the inverse size of the vortex core on the order of the local healing length of the condensate: 



2^2 



^ = 0.59 fi - -muj 2 {aRf . (399) 
The mean energy of O (384) is now a function of the position of the vortex core only, 

E = E no vortex + W{S) (400) 
where E no VO rtex is the energy of the condensate with no vortex and 



21n2 + l v(1q 2x 



(401) 



with v — 0.49312 and /i the chemical potential in the absence of vortex. This function 
W represents an effective potential seen by the vortex core. As shown in Fig. 14a this 
potential is maximal at the center of the trap so that it is actually an expelling potential 
for the vortex core: shifting the vortex core away from the center of the trap lowers the 
condensate energy. 

A method to stabilize the vortex is to rotate the harmonic trap around z at a fre- 
quency O (the trap is anisotropic in the x — y plane otherwise rotation would have 
no effect). Thermodynamical equilibrium will now be obtained in the frame rotating at 
the frequency O , where the harmonic trap is time independent. As this frame is non 
Galilean the Hamiltonian and therefore the Gross-Pitaevskii energy functional have to 
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be supplemented by the inert ial energy term — Q,L Z per atom, where L z is the angular 
momentum operator along z . The effective potential W(o) gets an extra term: 

W n (a) = W n=0 (3) - NMl(l - a 2 ) 2 (402) 

where Wq=q is the result (401) in the absence of rotation. As shown in Fig. 14b this extra 
term can trap the vortex core at the center of the harmonic trap if O is large enough. 

What happens if O is increased significantly ? It becomes favorable to put more 
vorticity in the condensate. As the vortices with charge larger than one are unstable the 
way out is to create several vortices with unit charge. This can be analyzed along the 
previous lines by a generalized multi- vortex ansatz, as discussed in [60]. A condensate 
with vortices has been recently obtained at the ENS in a rotating trap [61]. 

Another philosophy was followed at JILA: rather than relying on thermal equilibrium 
in a rotating trap to produce a vortex they used a "quantum engineering " technique 
[62] to directly induce the vortex by giving angular momentum to the atoms through 
coupling to electromagnetic fields [63] . It has also been suggested to imprint the phase of 
the vortex on the condensate through a lightshift induced by a laser beam whose spatial 
intensity profile has been conveniently tailored [64]. Such an imprinting technique has 
successfully led to the observation of dark and gray solitons in atomic condensates with 
repulsive interactions in Hannover [65] and in the group of W. Phillips at NIST. All these 
techniques illustrate again the powerfulness of atomic physics in its ability to manipulate 
a condensate. 

8 Phase coherence properties of Bose-Einstein con- 
densates 

Consider two Bose-Einstein condensates prepared in spatially well separated traps and 
that have 'never seen each other' (e.g. one rubidium condensate at JILA and one rubidium 
condensate at ENS). It is 'natural' to assume that these two condensates do not have a 
well defined relative phase. However the trend in the literature on Bose condensates is 
to assume that the two condensates are in a coherent state with a well defined relative 
phase, the so-called 'symmetry-breaking' point of view. So imagine that one lets the two 
condensates spatially overlap. Will interference fringes appear on the resulting atomic 
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Figure 14: In a 2D model, effective potential energy W of a vortex in a quasi-axisymmetric 
harmonic trap as function of the distance aR of the core from the trap center, for /jLq = SOHoj . 
(a) ft = and (b) fl = 0.045cl> . The unit of energy is NTiuj where uj is the oscillation 
frequency of the atoms in the trap. 



density or not ? 

One of the goals of this chapter is to answer this question and to reconcile the symmetry 
breaking point of view with the 'natural' point of view. 



8.1 Interference between two BECs 

At MIT a double well trapping potential was obtained by superimposing a sharp barrier 
induced with laser light on top of the usual harmonic trap produced with a magnetic 
field. In this way one can produce two Bose-Einstein condensates, one on each side of the 
barrier. The height of the barrier can be made much larger than the chemical potential 
of the gas so that coupling between the two condensates via tunneling through the wall 
is very small. In this way one can consider the two condensates as independent. 

One can then switch off the barrier and magnetic trap, let the two condensates bal- 
listically expand and spatially overlap. One then measures the spatial density of the 
cloud by absorption imaging. This spatial density exhibits clearly fringes [66] (see figure 
15). These fringes have to be interference fringes, as hydrodynamic effects (such as sound 
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waves) are excluded at the very low densities of the ballistically expanded condensates. 
We show here on a simple model that we indeed expect to see interference fringes in such 
an experiment, even if the two condensates have initially no well defined relative phase. 




Figure 15: Interference fringes between two condensates observed at MIT [66]. 



8.1.1 A very simple model 

In our simple modelization of an MIT-type interference experiment we will concentrate on 
the positions of the particles on an axis x connecting the two condensates so that we use 
a one-dimensional model enclosed in a box of size L with periodic boundary conditions. 
We assume that the system is initially in the Fock state 

l*> = if (403) 

with N/2 particles in the plane wave of momentum hk a and N/2 particles in the plane 
wave of momentum hkt : 

(x\k ajb ) = -T^exp [ik ajb x] . (404) 
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We assume that one detects the position of all the particles. What will be the outcome ? 
As the numbers of particles are exactly defined in the two modes a and b the relative 
phase between the atomic fields in the two modes is totally undefined. 

8.1.2 A trap to avoid 

If we calculate the mean density in the state given by (403) we find a uniform result 

(ft(x)$(x)) =N/L (405) 

and we may be tempted to conclude that no interference fringes will appear in the beating 
of two Fock state condensates. 

Actually this naive statement is wrong. Interference fringes appeared in a single 
realization of the experiment at MIT. We have therefore to consider the probability of the 
outcome of a particular density profile in a single realization of the measurement and not 
the average of the density profile over many realizations of the experiment. Indeed we 
will see that by interfering two independent Bose-Einstein condensates we get interference 
fringes on the density profile in each single realization of the experiment but the position 
of the interference pattern is random so that by averaging the density profile over many 
realizations we wash out the fringes. 

We wish to emphasize the following crucial point of the quantum theory: Whatever 
single-time measurement is performed on the system all the information about the out- 
comes of a single realization of the measurement procedure is contained in the N— body 
density matrix, here 

p=\q)(q\. (406) 

Indeed the only information we can get from quantum mechanics on a single realization 
outcome is its probability P , which can be obtained from p by 

P = Ti[6p] (407) 

where the operator O depends on the considered outcome. E.g. in our gedanken exper- 
iment P is the probability density of finding the N particles at positions xi, x 2 , ....x N 
and the operator O is expressed in terms of the field operator as 

6 = ^\x 1 )...^(x N )$(x N )..4(xi)- (408) 
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In a first quantized picture this corresponds to the fact that the probability density P is 
equal to the modulus squared of the N— body wavefunction. 

The complete calculation of the N— body distribution function P(xi, . . . ,x N ) for 
the state |\I>) in Eq.(403) is involved and we will see in the coming subsections how to 
circumvent the difficulty. But we can do a simple calculation of the pair distribution 
function of the atoms in state |^) : 



p(x u x 2 ) = {^\x l )^\x 2 )^{x 2 )i){x l )\^) 



(409) 
(410) 



We expand the field operator on the two modes 4> a:b and on other arbitrary orthogonal 
modes not relevant here as they are not populated in : 



$(x) = a(x\k a ) + b(x\k b ) + . . . 
where a and b annihilate a particle in state k a and k b respectively. We obtain 



(411) 



ip(x 2 )ip(xi)\^) = 



+ 



N (N 

N ( N 
2 V2 



N N 
(x 2 \k a )(xi\k a )\ — - 2 : k a ,— : k b ) 

i 1/2 N N 

(x 2 \k b )(xi\k b )\ — : k a , — - 2 : k b ) 



Nr i N N 

+ y [(X2\k a )(xi\k b ) + (x^ix^ka)^— - 1 : k a , — - 1 : A; 6 >.(412) 

The last line of this expression exhibits an interference effect between two amplitudes, 
that could not appear in the previous naive reasoning on the one-body density operator 
Eq.(405)! In the limit N ^> 1 and using the fact that the populated modes are plane 
waves the pair distribution function simplifies to 



p{x u x 2 ) ^ | 



1 + - cos 



[(ka - h)(xi - X 2 )] J ■ 



(413) 



This function exhibits oscillations around an average value equal to the square of the 
mean density. The oscillations are due to the interference effect in Eq.(412): they favor 
detections of pairs of particles with a distance \xi — x 2 \ equal to 2mi/\k a — k b \ (n 
integer) and they rarefy detections of pairs of particles with a distance (2n+l)7r/|/c a — k b \ . 
We therefore see on the pair distribution function a precursor of the interference fringes 
observed when the positions of all the particles are measured! 
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8.1.3 A Monte Carlo simulation 

By sampling the N— body distribution function P with a Monte Carlo technique, Ja- 
vanainen and Sung Mi Yoo in [67] made a numerical experiment with TV = 10 3 particles 
and k b — —k a . By distributing the measured positions in a given realization Xi, x 2 , ....x N 
among 30 position bins they obtained histograms like the ones in figure 16. It turns out 
that the density in the outcome of each realization of the numerical experiment can be 
fitted by a cosine: 

N 2 

where a and b are phases varying randomly from one realization to the other. In other 
words one has the impression that for each realization the system is in the state 



\0) N = 



|0) (415) 



with the angle = (0 a — b )/2 randomly distributed in [— 7r/2,7r/2] . Such a state, 
corresponding to a well defined phase between the two modes a and b , is called a phase 
state [68]. 

8.1.4 Analytical solution 

We wish to explain the result of the numerical experiment with an analytical argument. 
This has been done with slightly different points of view in [6, 69]. We give here what we 
think is the simplest possible presentation. 

Let us allow Poissonian fluctuations in the number of particles N a and N b , corre- 
sponding to the distribution probabilities: 



W) = ^e" w < e = a,b (416) 

with mean number of particles N a = N b = N/2 . These fluctuations become very small 
as compared to N when the number of particles becomes large: 

^ = -4= -> for N € oo. (417) 



The corresponding density operator is a statistical mixture of Fock states: 

oo 

p= Y, Pa{N a )Vt(N b )\N a :k a ,N b :k b ){N a :k a ,N b :k b \. (418) 

N a ,N b =0 
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Figure 16: For two different Monte Carlo realizations (a) and (b) of the gedanken experiment, 
histogram of the measured positions of N = 1000 particles for an initial Fock state with N/2 
particles in plane wave k a and N/2 particles in plane wave k^ = —k a [67]. The positions of 
the particles are expressed in units of 2n/(k a — k^) and are considered modulo 2n/(k a — . 

From this form one can imagine that a single realization of the experiment is in a Fock 
state, provided that one keeps in mind that N a and N b vary in an impredictable way 
from one experimental realization to the other. We known from the work [67] that there 
will be interference fringes in each experimental realization, but this fact is not intuitive. 

The same density operator can also be written in terms of a statistical mixture of 
phase states: 



From this form one can imagine that a single realization of the experiment is in a phase 
state, provided that one keeps in mind that the total number of particles N and the 
relative phase vary in an impredictable way from one realization to the other. This last 
form leads to the following algorithm to generate the positions of the particles according 
to the correct probability distribution: 

1. generate an integer N according to the Poisson distribution of parameter N 




(419) 
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p(x) = -!- e ikaX e ie + e ikbX e~ i6 2 . (421) 



2. generate according to a uniform probability distribution within — tt/2 and 7r/2 

3. generate the positions xi,....xn as if the system was in the state |0)jv , in which 
case all the particles are in the same single particle-state and the probability density 
P(xi, ....x N ) is factorized: 

N 

P(x u ....x N ) = l[p(x j ) (420) 
j=i 

where 

2L 

One then obtains interference fringes in each experimental realization, in a very explicit 
way. 

One could also use a third form of the same density operator p , that is a statistical 
mixture of Glauber coherent states: 

p t ^| C0h : iV^Woh : iW>)<coh : iV a 1/2 e*,coh : N^\. 

JO 27T JO 27T 

(422) 

This mathematical form is at the origin of the popular belief that condensates are in 
coherent states. From this form one can only imagine that a single realization of the 
experiment is in a coherent state, keeping in mind that the phases a and 6 b vary in an 
impredictable way from one realization to the other. In this representation the occurrence 
of interference fringes is straightforward. 

There is an important difference between the coherent states and the Fock or phase 
states: as the number of particles is a conserved quantity in the non-relativistic Hamil- 
tonian used to describe the experiments on atomic gases it seems difficult to produce 
a condensate in a coherent state in some mode tp , that is with p being a pure state 
|coh : a)(coh : a\ where a is a complex number. 

On the contrary one could imagine producing a condensate in a Fock state by mea- 
suring the number of particles in the condensate. One could then obtain a phase state by 
applying a 7r/2 Rabi pulse on the Fock state changing the internal atomic state a to a 
superposition (|a) + \ b))/y/2 where b is another atomic internal state; such a Rabi pulse 
has been demonstrated at JILA and has allowed the measurement of the coherence time 
of the relative phase between the a and b condensates [70]. 
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8.1.5 Moral of the story 

• there is in general no unique way of writing the density operator p as a statistical 
mixture. The canonical form corresponding to the diagonalization of p is always 
a possibility but not always the most convenient one. E.g. in our simple model the 
eigenbasis (Fock states) is less convenient than the non-orthogonal family of phase 
states (symmetry breaking states). 

• no measurement or no set of measurements performed on the system can distinguish 
between two different mathematical forms of the same density matrix as a statistical 
mixture. 

• the symmetry breaking point of view consists in writing (usually in an approximate 
way) the N— body density operator as a statistical mixture of Hartree-Fock states. 
One can then imagine that a given experimental realization of the system is a 
Hartree-Fock state, whose physical properties are immediate to understand as all 
the particles are in the same quantum state. 

• If the system is not in a state that is as simple as a Hartree-Fock state (e.g. in a Fock 
state for our simple model) it is dangerous to make reasonings on the single particle 
density operator (that is on the first order correlation function of the atomic field 
operator) to predict outcomes of single measurements on the system: the relevant 
information may be stored in higher order correlation functions of the field. 

8.2 What is the time evolution of an initial phase state ? 
8.2.1 Physical motivation 

Consider an interference experiment between two condensates A and B either in spa- 
tially separated traps or in different internal states (JILA-type configuration [70]). Assume 
that the two condensates have been prepared initially in a state with a well defined rel- 
ative phase ; this has actually been achieved at JILA. Let the system evolve freely for 
some time t . How long will the relative phase remain well defined ? This question is 
probably not an easy one to answer. We present here a simple model including only two 
modes of the field. In real life the other modes of the field are not negligible (see for 
example [71] for a discussion of finite temperature effects) and phenomena neglected here 
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such as losses of particles from the trap and fluctuations in the total number of particles 
may be important in a real experiment [7, 29]. 

We assume that the state of the system at time t = is a phase state. More 
specifically, expanding the N— th power in Eq.(415) with the binomial formula, we take 
as initial state: 

N ( N\ \ 1/2 

m = 0)) = 2""/ 2 E q [n^) j {Na - Nb) >a : 0a, N b : tf> b ) (423) 

where N b = N - N a and (j) a ^ b are the steady state condensate wavefunctions with N &ib 
particles in condensates A, B respectively. The time evolution during t is simple for 
each individual Fock states, as the system is then in a steady state with total energy 
E{N a ,N b ): 

\N a : a , N b : fa) -+ e ~ iE ^ N ^ n \N a : a , N b : fa). (424) 

The time evolution of the phase state Eq.(423) is much more complicated: the state vector 
\^(t)) is a sum of many oscillating functions of time. 

8.2.2 A quadratic approximation for the energy 

The discussion can be greatly simplified if one uses the fact that the binomial factor in 
Eq.(423) for large N is a function of N a and N b sharply peaked around N a = N b = N/2 
with a width y/N : from Stirling's formula n\ ~ (n/e) n y/27rn we obtain indeed 

1 Nl 1 r N ^ e - Na \o g (N a /N)-N b log(N b /N) _ J_V /2 g -(iV a -iV 6 ) 2 /(2iV) 



2 N N a \N b \ ^/2tt2 n \N a N b J UN, 

(425) 

We therefore expand the energy E in powers of N a - N/2 and N b - N/2 up to second 
order. 

E(N a ,N b ) = E{N/2,N/2) + (N a -N/2)d Na E+(N b -N/2)d Nh E 



+\{N a - N/2fdlE + \{N b - N/2fdl b E 
+ {N a - N/2)(N b - N/2)d Na d Nb E + . . . , (426) 



all the derivatives being taken in (N a , N b ) = (N/2, N/2) . Note that the first derivatives 
of the energy are the chemical potentials fi ajb of the two condensates; as the condensates 
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are independent condensates (there is no mechanism locking the relative phase of the 
condensates) one has in general /i a ^ \i b . As we restrict to the set of occupation numbers 
such that N a + N b = N we can rewrite the expansion of the energy using N a - N/2 = 
-(N b -N/2) = (N a -N b )/2: 

E(N a ,N b ) ~ E(N/2,N/2) + ^(A a - N b )(fi a - + ^(N a - N b ) 2 X (427) 
where we have introduced the quantity 

x = -1 \(d N -d Nh ) 2 E\ . (428) 

2k L b lN n =N h =N/2 K } 



8.2.3 State vector at time t 



If one uses the quadratic approximation of the energy the system evolves from the initial 
state Eq.(423) to the state 

JV / iu| \ V 2 

|*(t)> = £ -f— ^We-W-W*^ : ^ . ^ (42 9) 

The contribution of the term linear in iV a — N b in Eq.(427) is contained in the quantity 

v=^ b -f^ a ). (430) 

The resulting effect on the time evolution is simply to shift the relative phase between the 
condensate from 9 to 9 + vt : this is a mere phase drift with a velocity v . This phase 
drift takes place only if the 'frequencies' fi a /h and fi b /h of the atomic fields in A and 
in B are different. 

The effect of the quadratic term in Eq.(427) is to spread the relative phase of the 
two condensates. This effect is formalized in [6], we give here the intuitive result. The 
spreading of the phase can be understood in analogy with the spreading of the wavepacket 
of a fictitious massive particle, with the relative phase 9 being the position x of the 
particle and the occupation number difference — N a being the wavevector k of the 
particle. The energy term proportional to % plays the role of the kinetic energy of the 
particle responsible for the spreading in position. The effective mass of the fictitious 
particle is M such that 

^Na-N.fx^^ (431) 
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so that 

M = ^. (432) 

M 

Replacing the discrete sum in Eq.(429) by an integral we formally obtain the expansion 
of the time dependent state vector of the fictitious particle over the plane waves in free 
space. In this case the variance of the position of the fictitious particle spreads as 

Ax 2 W = Ax 2 (0)+(^)V (433) 

Within the approximation (425) the wavepacket of the fictitious particle is a Gaussian in 
momentum space, with a standard deviation Ak = (N/2) 1 / 2 . Initially the position x is 
well defined with a spread ~ 1/Ak <C 1 . The relative phase of the condensates will start 
becoming undefined when the position spread Ax of the fictitious particle becomes on 
the order of unity. This happens after a time 

M 2y/2 

Vead ~ ^ = (434) 

At times much longer than £ spre ad it is not correct to replace the discrete sum over 
(N a - N b )/2 by an integral. The discreteness of N a - N b leads to reconstructions of a 
phase state (the so-called revivals) at times t q = qn/x , Q integer: one can check indeed 
from Eq.(429) that a phase state is reconstructed with a relative phase O + vtq + qir /2 for 
N even and + vt q for N odd. The observability of even the first revival at time ti is 
a non trivial question: the revivals are easily destroyed by decoherence phenomena such 
as the loss of a few particles out of the condensate due to inelastic atomic collisions [7], 
and effects of the non-condensed fraction also need to be investigated. This fragility of 
the revivals is not surprising if one realizes that the state vector \^(t)) in Eq.(429) is a 
Schrodinger cat at time ti/2 , that is a coherent superposition of the N particles in some 
state 0i and of the N particles in some state </5> 2 orthogonal to <j>i : the revival at time 
ti is suppressed if the Schrodinger cat at time ti/2 is transformed by decoherence into 
a statistical mixture of the states \N : 0i ;2 ) , which is difficult to avoid for large values of 
N (see the lecture notes of Michel Brune in this volume). 

8.2.4 An indicator of phase coherence 

To characterize the degree of phase correlation between the two condensates it is natural 
to consider the average of (cftb) where a, b annihilate a particle in condensates A and 
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B respectively. Consider indeed the average over many experimental realizations of some 
one-body observable O sensitive to the relative phase of the two condensates. This 
observable necessarily has a non- vanishing matrix element between the modes (j) a and 
4> b so that in second quantized form the part of (O) sensitive to the relative phase involves 
(atb) . E.g. in the case of spatially separated condensates one can beat on a 50 - 50 matter 
waves beam splitter atoms leaking out of the condensates and detect the atoms in the 
output channels of the beam splitter [6] ; the number of counts in the + output channel 
averaged of many experimental realizations is proportional to the expectation value of 

Expanding this product of operators we get 'diagonal' terms such as a^a not sensitive to 
the relative phase, and crossed terms (actually interference terms!) such as cftb sensitive 
to the phase. In the JILA-type configuration, where the condensates A and B are in 
different internal atomic states, an observable O similar to Eq.(435) has been achieved 
by mixing the internal states of the two condensates by a 7r/2 electromagnetic pulse and 
by measuring the mean density of atoms in A and B [70] . 

From the Schwartz inequality \{u\v)\ < \\u\ \ \ \v\\ and setting \u) = a\^) , \v) = b\^) 
we obtain an upper bound for the expectation value of atb : 

< (^la^l^) 1 / 2 ^!^^) 1 / 2 . (436) 

The case of a maximally well defined relative phase corresponds to an equality in this 
inequality, obtained if \u) and \v) are proportional. In the present situation of equal 
mean numbers of particles N/2 in A and in B this corresponds to being a phase 
state. 

For an initial phase state it is possible to calculate the expectation value of ttfb as 
function of time from the expansion (429). One obtains after simple transformations the 
sum 

(a r b)(t) = — e -^ 6 + vt ) V ( N e ixt[N a -(N b -i)] (m) 

with N b = N - N a as in Eq.(429). After inspection one realizes that this sum is the 
binomial expansion of a (N - 1)- th power so that the final result is [72]: 

(atg) (t) = ^ e - 2 ^) cos^- 1 X t. (438) 
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From this very simple expression one can calculate the time t c after which the relative 
phase has experienced a significant spread. For short times %t <C 1 one can expand the 
cosine function in Eq.(438) to second order in t : 

cos^ X t = e Nl0 % C0Sxt ~ e"^**) 2 / 2 . (439) 

One obtains a Gaussian decay of phase coherence with a collapse time 

te = W^ (440) 

equivalent to the rougher estimate Eq. (434) up to a numerical factor. One can also easily 
see the revivals (reconstruction of |\I>) to a phase state) at times t q = qn/x when the 
cosine function is equal to ±1 in Eq.(438). 

Formula (440) can be used to calculate the coherence time of the relative phase of the 
condensates in the present zero-temperature model. As an interesting application of this 
formula we now show that the spreading time of the relative phase can be significantly 
different for mutually interacting and non- mutually interacting condensates. Assume for 
simplicity that the two condensates are stored in cubic boxes of identical size L and with 
periodic boundary conditions. In the MIT-type configuration the two boxes are spatially 
separated and the atoms are in the same internal state; the energy of a configuration with 
N a ,N b atoms in the condensates A,B is then 



E = 



2L 3 L 



N 2 a +Nf . (441) 



From Eq.(428) this form of E leads for an initial phase state to a collapse time of the 
relative phase 

t c = N^ 2 — (442) 

where p = N/ (2L 3 ) is the mean spatial density in each of the condensates. In the JILA- 
type configuration the atoms are in the same spatial box but in different internal states; 
the energy of a configuration with N a ,N b atoms in the condensates A, B is given now 
by Eq.(294) if the two internal states are subject to spatial demixing, or by Eq.(295) if 
there is no demixing instability. The collapse time is then given by 

*c = ^ l/2 ~ r ^7 \T7oi (443) 

P[9aa + gbb -^{9aa9bb) 1/2 \ 
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for a demixed condensates and by 

t c = iV 1 / 2 - — (444) 

P[9aa + 9bb ~ 2gab\ 

for fully overlapping condensates. When the coupling constants among the various internal 
states are close to each other the denominator in Eqs. (443,444) can become small, which 
results in a relative phase coherence time t c much larger than in the MIT-configuration 
Eq.(442). This fortunate feature of close coupling constants is present for rubidium in the 
JILA experiment [70]! 

In real life the condensates are usually stored in harmonic traps; the simple formulas 
obtained for a cubic box have to be revisited. This has been done analytically for spatially 
separated condensates [73, 74] and numerically for mutually interacting condensates [29, 
75]. 



9 Symmetry breaking description of condensates 

We have already seen in chapter 8 that it is very convenient, physically, to introduce 
phase states to understand the phenomenon of interference between two Bose-Einstein 
condensates: rather than assuming that two Bose-Einstein condensates that "have never 
seen each other" are in Fock states, one assumes that they are in a phase state with a 
relative phase varying in an unpredictable way for any new experimental realization. One 
can even suppose that the condensates are in coherent states of the atomic field; this 
description is said to 'break the symmetry', here the U(l) symmetry associated to the 
invariance of the Hamiltonian by a change of the phase of the atomic field operator. 

In this chapter we consider other examples of symmetry breaking descriptions: SO (3) 
symmetry breaking (case of spinor condensates) and spatial translational symmetry break- 
ing (case of one dimensional condensates with attractive interactions). In both cases the 
procedure is the same: the ground state of the system is symmetric, its mean-field ap- 
proximation by Hartree-Fock states breaks the symmetry. In both cases we will consider 
Gedanken experiments whose single outcomes can be predicted from the exact ground 
state and from the Hartree-Fock state. This will illustrate the ability of the mean-field 
approximation to allow physical predictions in an easy and transparent way, correct in 
the limit of a large number of particles. 
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9.1 The ground state of spinor condensates 

The alkali atoms used in the Bose- Einstein condensation experiments have an hyperfine 
structure in the ground state, each hyperfine level having several Zeeman sublevels. We 
have up to now ignored this structure in the lecture, as we were implicitly assuming that 
the atoms were polarized in a well defined Zeeman sublevel. 

Consider for example 23 Na atoms used at MIT in the group of Wolfgang Ketterle. 
The ground state has an hyperfine splitting between the lower multiplicity of angular 
momentum F = 1 and the higher multiplicity of angular momentum F = 2 . All the 
three Zeeman sublevels m F = 0, ±1 of the lower multiplicity F = 1 cannot be trapped 
in a magnetic trap (if mp — —1 is trapped than mp — +1 which experiences an 
opposite Zeeman shift is antitrapped) . But they can all be trapped in an optical dipole 
trap, produced with a far off-resonance laser beam, as the Zeeman sublevels experience 
then all the same lightshift. This optical trapping was performed at MIT [76], opening 
the way to a series of interesting experiments with condensates of particles of spin one 
[77]. 

We concentrate here on a specific aspect, the ground state of the spinor condensate, 
assuming for simplicity that the atoms are stored in a cubic box with periodic boundary 
conditions. 

9.1.1 A model interaction potential 

We have to generalize the model scalar pseudo-potential of Eq.(73) to the case of particles 
having a spin different from zero. As we want to keep the simplicity of a contact interaction 
potential we choose the simple form 

y(l,2)=V spin (l,2)^(rl-r1) 

that is the product of an operator acting only on the spin of the particles 1 and 2, and 
of the usual regularized contact interaction acting only on the relative motion of the two 
particles. The interaction potential V(l,2) has to be invariant by a simultaneous rotation 
of the spin variables and of the position variables of the two particles. As the contact 
interaction is already rotationally invariant, the spin part of the interaction V sp i n (l,2) 
has to be invariant by any simultaneous rotation of the two spins. 



d 



dr 



12 



(ri2 ■ ) 



(445) 
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This condition of rotational invariance of V sp i n (l, 2) is easy to express in the coupled 
basis obtained by the addition of the two spins of particle 1 and particle 2: within each 
subspace of well defined total angular momentum V sp i n (l,2) has to be a scalar. Let us 
restrict to the case studied at MIT, with spin one particles. By addition of F = 1 and 
F — 1 we obtain a total angular momentum F tot = 2 , 1 or 0, so that one can write 

V spin (l, 2) = g 2 P Ftot=2 (l, 2) + 5iPF tot =i(l, 2) + 0,^=0(1, 2) (446) 

where the g 's are coupling constants and the F(l,2) 's are projectors on the subspace 
of particles 1 and 2 with a well defined total angular momentum F tot . At this stage we 
can play a little trick, using the fact that the states of F to t = 1 are antisymmetric by the 
exchange of particles 1 and 2 (whereas the other subspaces are symmetric). The regu- 
larized contact interaction scatters only in the 5 -wave, where the external wavefunction 
of atoms 1 and 2 is even by the exchange of the positions fi and r 2 ; as our atoms are 
bosons, the spin part has also to be symmetric by exchange of the spins of atoms 1 and 2 
so that the 'fermionic' part of V sp i n (l, 2) , that is in the subspace F tot = 1 , has no effect. 
We can therefore change gi at will without affecting the interactions between bosons. 
The most convenient choice is to set g x = g 2 so that we obtain 

V spin (l, 2) = <feld(l, 2) + {g Q - 92)P Ftot =o(h 2) (447) 

where Id is the identity. The subspace F tot = is actually of dimension one, and it is 
spanned by the vanishing total angular momentum state \ipo(l, 2)) . Using the standard 
basis \m = —1, 0, +1) of single particle angular momentum with z as quantization axis, 
one can write 

|Vo(l, 2)> = ~ [\ + 1, -1> + | - 1, +1) - |0, 0)] . (448) 

A more symmetric writing is obtained in the single particle angular momentum basis 
\xjjjjz) used in chemistry, defined by 

| + 1) = -±(\ X )+i\y)) (449) 

|-1> = +^(\x)-i\v)) (450) 
|0) = \z). (451) 
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The vector \a) in this basis (a = x,y, z) is an eigenvector of angular momentum along 
axis a with the eigenvalue zero. One then obtains 

\Mh 2)) = [\x, x) + \y, y) + \z, z)\ . (452) 

To summarize the part of the Hamiltonian describing the interactions between the 
particles can be written, if one forgets for simplicity the regularizing operator in the 
pseudo-potential: 



TJ 92 
ttmt = — 



a,P=x,y,z 

Jd 3 r £ faftMp. (453) 



2 , R 

a,0=x,y,z 

,90-92 



6 

a,p=x,y,z 



where ij> a (r) is the atomic field operator for the spin state \a) . This model Hamiltonian 
has also been proposed by [78, 79, 80]. 

9.1.2 Ground state in the Hartree-Fock approximation 

As we are mainly interested in the spin contribution to the energy we assume for simplicity 
that the condensate is in a cubic box of size L with periodic boundary conditions. We 
assume that the interactions between the atoms are repulsive (g2,9o > ) and we suppose 
that there is no magnetic field applied to the sample. 

We now minimize the energy of the condensate within the Hartree-Fock trial stat- 
evectors |A^ : <£) with the constraint that the number of particles A^ is fixed ( is 
normalized to unity) but without any constraint on the total angular momentum of the 
spins. The external part of the condensate wavefunction is simply the plane wave with 
momentum k = whereas the spinor part of the wavefunction remains to be determined: 

W> = 75/2 E C M with El c a | 2 = 1 - (454) 

a=x,y,z a 

Using the model interaction Hamiltonian Eq.(453) we find for the mean energy per particle 
in the condensate 

E Nq-1 N -l, ,, rr , 

N-r^ 92+ ^ i9 °- g2M (455) 
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where we have introduced the complex quantity 

A= £ 4 = c 2 (456) 

a=x,y,z 

where c is the vector of components (c x , c y , c z ) . We have to minimize the mean energy 
over the state of the spinor. 

• Case g 2 > go 

This is the case of sodium [77]. As the coefficient g — g 2 is negative in Eq.(455) we have 
to maximize the modulus of the complex quantity A . As the modulus of a sum is less 
than the sum of the moduli we immediately get the upper bound 

\A\< £ |c | 2 = l (457) 

a=x,y,z 

leading to the minimal energy per particle 

E N -I , N -l , 

The upper bound for \A\ is reached only if all complex numbers c 2 a have the same phase 
modulo 27T . This means that one can write 

c a = e i6 n a (459) 

where 6 is a constant phase and n=(n x ,n y ,n z ) is any unit vector with real components. 
Physically this corresponds to a spinor condensate wavefunction being the zero angular 
momentum state for a quantization axis pointing in the direction of ft . The direction 
of n is well defined in the Hartree-Fock ansatz, but it is arbitrary as no spin direction 
is privileged by the Hamiltonian. We are facing symmetry breaking, here a rotational 
SO (3) symmetry breaking, as we shall see. 

• Case g 2 < go 

In this case we have to minimize \A\ to get the minimum of energy. The minimal value 
of \A\ is simply zero, corresponding to spin configurations such that 

c 2 = £ cl = (460) 

a=x,y,z 
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with an energy per condensate particle 

N = -W^ 9 *- (461) 
To get more physical understanding we split the vector c as 

?=R + iT (462) 

where the vectors R and / have purely real components. Expressing the fact that 
the real part and imaginary part of c 2 vanish, and using the normalization condition 
c • c* = 1 in Eq.(454) we finally obtain 

R-T = (463) 
R 2 = P = \. (464) 

This means that the complex vector c is circularly polarized with respect to the axis Z 
orthogonal to / and R . Physically this corresponds to a spinor condensate wavefunction 
having an angular momentum ±h along the axis Z . The direction of axis Z is well 
defined in the Hartree-Fock ansatz but it is arbitrary. 



9.1.3 Exact ground state of the spinor part of the problem 

Imagine that we perform some intermediate approximation, assuming that the particles 
are all in the ground state k — of the box but not assuming that they are all in the 
same spin state. We then have to diagonalize the model Hamiltonian 

tfspin = a^apaa + 773 (po - 92) A* A (465) 

a,0=x,y,z 

where a a annihilates a particle in state \k = 0)\a) (a = x,y, z) and where we have 
introduced 

A = &l + al + al (466) 

Up to a numerical factor A annihilates a pair of particles in the two-particle spin state 
I -00 ( 1 5 2) ) of vanishing total angular momentum, as shown by Eq.(452). 

The Hamiltonian Eq.(465) can be diagonalized exactly [81]. This is not surprising as 
(i) it is rotationally invariant and (ii) the bosonic A^— particle states with a well defined 
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total angular momentum S No can be calculated: one finds that S No — N Q , N — 2, . . . , 
leading to degenerate multiplicities of H spin of degeneracy 25^ + 1 . 

In practice one may use the following tricks: The double sum proportional to g 2 in 
Eq.(465) can be expressed in terms of the operator number N Q of condensate particles 
only, 

a 

So diagonalizing H spin amounts to diagonalizing A^A ! 

Second the total momentum operator S of the N spins, defined as the sum of all 
the spin operators of the individual atoms in units of % , can be checked to satisfy the 
identity 

S • S + i f i = N (N + 1) (468) 
so that the Hamiltonian for N Q particles becomes a function of S [81]: 

ffspin = j^MNo - 1) + ^(9o - 92) [MNv + 1) - I • I] . (469) 

We recall that S • S = Sn (Sn + 1) within the subspace of total spin Sn ■ 

When g 2 < go the ground state of i^ spin corresponds to the multiplicity S No = , 
containing e.g. the state with all the spins in the state |+) . In this case the A^ — particle 
states obtained with the Hartree-Fock approximation are exact eigenstates of H spin . 

When g 2 > go the ground state of H sp - m corresponds to the multiplicity of minimal 
total angular momentum, S No — 1 for odd or S No — for even. In this case 
the Hartree-Fock state is a symmetry breaking approximation of the exact ground state 
of -ffspin ■ The error on the energy per particle tends to zero in the thermodynamical 
limit; for N even one finds indeed 

f = -3Z*(*-*>- ( 47 °) 

But what happens if one restores the broken symmetry by summing up the Hartree- 
Fock ansatz over the direction ft defined in Eq.(459)? Assume that N is even; one 
has then to reconstruct from the Hartree-Fock ansatz a rotationally invariant state. This 
amounts to considering the following normalized state for the AT spins: 

\V) = jN + lj—\N :fi) (471) 
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where dPn indicates the integration over the unit sphere (that is over all solid angles) 
and \N Q : ft) is the state with N Q particles in the single particle state 

\n) = n x \x) + n y \y) + n z \z). (472) 

The state vector , being non zero and having a vanishing total angular momentum, 
is equal to the exact ground state of H spin ! 

The expression (471) can be used as a starting point to obtain various forms of |\I>) . 
If one expresses the Hartree-Fock state as the N -th power of the creation operator 
Ed ®a n a acting on the vacuum |vac) , and if one expands this power with the usual 
binomial formula, the integral over ft can be calculated explicitly term by term and one 
obtains: 

I*) =Af(A*) No/2 | vac) (473) 

where Af is a normalization factor and the operator A is defined in Eq.(466). Formula 
(473) indicates that |\I>) is simply a 'condensate' of pairs in the pair state |V>o(l, 2)) . It 
can be used to expand |\I>) over Fock states with a well defined number of particles in 
the modes m = 0,m = ±1 , reproducing Eq.(13) of [81]. 

To be complete we mention another way of constructing the exact eigenvectors and 
energy spectrum of H spin . The idea is to diagonalize A^A using the fact that A obeys 
a commutation relation that is reminiscent of that of an annihilation operator: 

[i,it] = 47V + 6. (474) 

In this way A^ acts as a raising operator: acting on an eigenstate of A^A with eigenvalue 
A and N particles, it gives an eigenstate of A^A with eigenvalue \-\-4N + 6 and with 
Nq + 2 particles. One can also check from the identity (468) that the action of A^ does 
not change the total spin: 

[A\S-S}= 0. (475) 

By repeated actions of A^ starting from the vacuum one arrives at Eq.(473), creating 
the eigenstates with AT even and vanishing total spin S — . By repeated actions of 
A^ starting from the eigenstates with N = 2 and total spin S = 2 (e.g. the state 
| + +) ) one obtains all the states with A^ even and total spin S = 2 . More generally 
the eigenstate of H sp i n with total spin S , a spin component m = S along z and AT 
particles is: 

\\N Q ,S,m = S) cx ^) {N °~ S)/2 \S : +1) (476) 
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where \S : +1) represents S particles in the state | + 1) . From Eq.(476) one can 
generate the states with spin components m — S — 1, . . . , —S by repeated actions of the 
spin-lowering operator S- = S x — iS y in the usual way. We note that formula (476) was 
derived independently in [82]. 

9.1.4 Advantage of a symmetry breaking description 

Imagine that we have prepared a condensate of sodium atoms ( g 2 > go ) in the collective 
ground spin state, and that we let the atoms leak one by one out of the trap, in a way 
that does not perturb their spin. We then measure the spin component along z of the 
outgoing atoms. Suppose that we have performed this measurement on k atoms, with 
k <C AT ■ We then raise the simple question: what is the probability pk that all the k 
detections give a vanishing angular momentum along z ? 

Let us start with a naive reasoning based on the one-body density matrix of the 
condensate (even if the reader has been warned already in §8.1.2 on the dangers of such 
an approach!). The mean occupation numbers of the single particle spin states \m = — 1) , 
\m = 0) and \m = +1) in the initial condensate are obviously all equal to N Q /3 , as the 
condensate is initially in a rotationally symmetric state. The probability of detecting the 
first leaking atom in \m = 0) is therefore 1/3 . Naively we assume that since k < N Q 
the detections have a very weak effect on the state of the condensate and the probability 
of detecting the n -th atom ( n < k ) in the m = channel is nearly independent of the 
n-1 previous detection results. The probability for k detections in the m = channel 
should then be 

P?™ = ^- (477) 
Actually this naive reasoning is wrong (and by far) as soon as k > 2 . The first 
detection of an atom in the m — channel projects the spin state of the remaining 
atoms in 

|* 1 >=A/iao|*> (478) 
where a Q annihilates an atom in spin state m = , |^) is the collective spin ground 
state (471) and M is a normalization factor. The probability of detecting the second 
atom in m = (knowing that the first atom was detected in m = ) is then given by 

P2 (^il^ol^i) 



Pi (*l|Em=-l^m|*l)' 



(479) 
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The denominator is simply equal to N Q — 1 as |\I>i) is a state with N Q — 1 particles. Using 
the integral form (471) and the simple effect of an annihilation operator on a Hartree-Fock 
state, e.g. 

a 2 \N : ft) = [N (N - 1)] 1/2 ^|AT - 2 : ft) (480) 
we are able to express the probability in terms of integrals over solid angles: 

/ d 2 n [ d 2 n' n 2 z n'*{n.n') N «- 2 
Pi - I 1 111 (481) 



Pi J d 2 n J tfn' n z n' z (n-n') 



t\N -l 



We suggest the following procedure to calculate these integrals. One first integrates over 
ft' for a fixed ft , using spherical coordinates relative to the 'vertical' axis directed along 
ft: the polar angle 9' is then the angle between ft' and ft so that one has simply 
ft -ft' = cos0' . The integral over 9' and over the azimuthal angle <f>' can be performed, 
giving a result involving only n z . The remaining integral over ft is performed with the 
spherical coordinates of vertical axis z . This leads to 

H + ^rry <«> 

The ratio P2/P1 is therefore different from the naive (and wrong!) prediction (477). 
For No = 2 one finds P2/P1 = 1 so that the second atom is surely in m = if the first 
atom was detected in m = . As the two atoms were initially in the state with total 
angular momentum zero, this result could be expected from the expression (448) of the 
two-particle spin state. In the limit of large A^ we find that once the first atom has been 
detected in the m — channel, the probability for detecting the second atom in the same 
channel m = is 3/5 . This somehow counter-intuitive result shows that the successive 
detection probabilities are strongly correlated in the case of the spin state (471). 

The exact calculation of the ratio 

/ d 2 ft [ d 2 n' n k z +1 n' z k+1 (n • ')*>-(*+!) 
Pk^ri _ J J ( 483 ) 



Pk I d 2 n J d 2 n' n k z n' z k (n • ff')^ 0- * 



is getting more difficult when k increases. The large limit for a fixed k is easier to 
obtain: in the integral over ft' the function (n • n') N °~( k+1 ^ is extremely peaked around 
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ft' — ft so that we can replace nf +l by n k+1 . This leads to 

lim £*±i = ^±l. (484) 

We now give the reasoning in the symmetry breaking point of view, which assumes 
that a single experimental realization of the condensate corresponds to a Hartree-Fock 
state \N : n) with the direction ft being an impredictable random variable with uniform 
distribution over the sphere. If the system is initially in the spin state |iVo : ft) there 
is no correlation between the spins, and the probability of having k detections in the 
channel m = is simply (n 2 z ) h . One has to average over the unknown direction ft to 
obtain 

pf = f d ^ n f = -i-. (485) 
Fk J 4tt z 2& + 1 v } 

One recovers in an easy calculation the large iVo limit of the exact result, Eq.(484)! We 

note that the result (485) is much larger than the naive (and wrong) result (477) as soon 

as k ^> 1 . 



9.2 Solitonic condensates 

We consider in this section a Bose-Einstein condensate with effective attractive inter- 
actions subject to a strong confinement in the x — y plane so that it constitutes an 
approximate one-dimensional interacting Bose gas along z . Such a situation is interest- 
ing physically as it gives rise in free space to the formation of 'bright' solitons well known 
in optics but not yet observed with atoms. Also the model of a one-dimensional Bose gas 
with a S interaction potential has known exact solutions in free space, that can be used 
to test the translational symmetry breaking Hartree-Fock approximation. 



9.2.1 How to make a solitonic condensate ? 

Consider a steady state condensate with effective attractive interactions in a three di- 
mensional harmonic trap. The confinement in the x — y plane is such that the trans- 
verse quanta of oscillation %uj x , y are much larger than the typical mean field energy per 
particle 7Vo|^||0| 2 , where is the condensate wavefunction with particles. This 
confinement prevents the occurrence of a spatial collapse of the condensate (see §5.2.1). 
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The confinement is however not strong enough to violate the validity condition of the 
Born approximation for the pseudo-potential, k\a\ <C 1 with k ~ (mojx^/h) 1 / 2 . 

In this case we face a quasi one-dimensional situation, where the condensate wave- 
function is approximately factorized as 



(/>(x,y,z) = ip(z)xx{x)x y (y) 



(486) 



where Xx and Xy are the normalized ground states of the harmonic oscillator along x 
and along y respectively. By inserting the factorized form (486) in the Gross- Pit aevskii 
energy functional Eq.(139) and by integrating over the directions x and y we obtain an 
energy functional for ip : 



E[il>,fl>*]=N Q J dz 



~h 2 


dip 


2m 


dz 



+ 1 -mu } y\4,(z)\ 2 + l -N o9l M{z)\ i 



(487) 



where we have dropped the zero-point energy of the transverse motion and we have called 
gid the quantity 



9u = gJdxJ dy MaOlUWI 4 = O ^^f' 2 - 



The corresponding time independent Gross-Pitaevskii equation for ip is 



^mu 2 z 2 + N g ld \<iP(z)\ 2 i/>(z). 



(488) 



(489) 



The energy functional Eq.(487) corresponds to a one-dimensional interacting Bose gas 
with an effective coupling constant between the atoms equal to gu , that is one can 
imagine that the particles have a binary contact interaction 



V(z u z 2 ) = gidd(zi ~ z 2 ). 



(490) 



Note that such a Dirac interaction potential leads to a perfectly well defined scattering 
problem in one dimension, contrarily to the three dimensional case. 

Imagine now that we slowly decrease the trap frequency along z while keeping intact 
the transverse trap frequencies, until uj z vanishes. What will happen then? If g was pos- 
itive the cloud would simply expand without limit along z . With attractive interaction 
the situation is dramatically different: due to the slow evolution of uj z the condensate 
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wavefunction will follow adiabatically the minimal energy solution of the Gross- Pit aevskii 
equation. For u z — this minimal energy solution is the so-called bright soliton, well 
known in non-linear optics. We recall the analytic form of the solitonic wavefunction: 

= ^JT) ^ 

where I is the spatial radius of the soliton: 

i = --f?— . (492) 
N mg ld 

Note that this size I results of a compromise between minimization of kinetic energy by 
an increase of the size and minimization of interaction energy by a decrease of the size, so 
that the typical kinetic energy per particle h 2 / (ml 2 ) is roughly opposite to the interaction 
energy per particle N gi d /l . We also give the corresponding chemical potential: 



M = -^ 2 ^. (493) 



We briefly address the validity of the Gross- Pit aevskii solution (491). As we have 
pointed out in the three dimensional case (see for example §3.2.1) we wish that the Born 
approximation for the interaction potential be valid. In one dimension the S interaction 
potential can be treated in the Born approximation only if the relative wave vector of the 
colliding particles is high enough (in contrast to the three-dimensional case) : 



h 2 k 



mgid 



> 1. (494) 



This condition can be obtained of course from a direct calculation, but also from a dimen- 
sionality argument ( mgu/h 2 is the inverse of a length) and from the fact that the Born 
approximation should apply in the limit gia — > for a fixed k . If we use the estimate 
k ~ l/l we obtain the condition 

~ No > 1, (495) 



mgi d l 

implicitly valid here as we started from a condensate! 

Another phenomenon neglected in the prediction (491) is the spreading of the center 
of mass coordinate during the switch-off of the trapping potential along z . Whereas 
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Eq.(491) assumes that the abscissa of the center of the soliton z Q is exactly the 
spreading of the center of mass leads in real life to a finite width probability distribution 
for zq . This spreading can be calculated simply for an almost pure condensate N ~ N , 
using the fact that the center of mass coordinate operator Z and the total momentum 
operator P of the gas along z axis are decoupled from the relative coordinates of the 
particles in a harmonic potential, in presence of interactions depending only on the relative 
coordinates. To prove this assertion one expresses the operators Z and P in terms of 
the position and momentum operators of each particle i of the gas: 

Z = ^X> (496) 

JV i=l 

P = X> (497) 

and one derives the following equations of motion in Heisenberg point of view: 

dZ _ P 

dt ~ Nm 
dP 

a — = -NmuZ(fi)Z. (499) 

The spreading acquired by Z is not negligible when it becomes comparable to the size I 
of the soliton. 

The spreading of Z is interesting to calculate in the absence of harmonic confinement 
along z , uj z = , with the simple assumption that all the particles of the gas are at time 
t = in the soliton state of Eq.(491). As P is a constant of motion for uj z = 
one has simply 

m = Z(0) + ^ (500) 

so that the variance of the center of mass coordinate at time t is 

Var(Z)(t) = Var(Z)(0) + ^<^(«)^ + PZ(0)) + ^Vax(P). (501) 

One then replaces Z(0) and P by the sums (496, 497). As the single particle wave- 
function ip has vanishing mean position and mean momentum all the 'crossed terms' 



(498) 
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expectation values involving two different particles vanish. As ip(z) is a real wavefunc- 
tion one finds also (ip\zp -\- pz\ip) = so that the contribution linear in time vanishes. 
One is left with 

Var(Z)(t) = + ^<WlV->- ( 502 ) 

The variance of Z , initially TV times smaller than the single particle variance (ip\z 2 \tp) , 
becomes equal to the single particle variance after a time 

where we used the explicit expressions 

Wz 2 W = ^ (504) 
WfW = |r- (505) 

The spreading phenomenon of the position of the soliton is formally equivalent to the 
spreading of the relative phase of two condensates initially prepared in a phase state (see 
§8.2). The critical time t c in (503) scales as N^ 2 h/\fi\ as in Eq.(442). 

9.2.2 Ground state of the one- dimensional attractive Bose gas 

We consider here the model of the one-dimensional gas of N bosonic particles interacting 
with the contact potential Eq.(490) and in the absence of any confining potential. 

It turns out that in this model with gu > one can calculate exactly the eigenenergies 
and eigenstates of the Hamiltonian for N particles using the Bethe ansatz [83]. We 
consider here the less studied attractive case gi d < , where several exact results are also 
available. In particular the exact expression for the ground state energy is known [84]: 



1 m 9id 



E Q (N) = --^N(N*- 1) (506) 
and the corresponding N— particle wavefunction of the ground state is [85]: 



W(zi,...,z N ) =A/"exp 



mgi d „ 

9 t2 2^ \ Z i Z 3 

l<i<j<N 



(507) 
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To determine the normalization factor Af we enclose the gas in a fictitious box of size L 
tending to +00 : 5 

To what extent can we recover these results using a Hartree-Fock ansatz \N : ip) for 
the ground state wavefunction? As discussed around Eq.(151) we get a mean energy for 
the Hartree-Fock state very similar to Eq.(487): 



E[ip,il>*] = N Jdz 



~h 2 


dip 


2m 


dz 



+ -{N-l)g ld \^{z)\ l 



(509) 



We minimize this functional using the results of §9.2.1, replacing N Q by N — 1 , and we 
obtain 

1 Wxd 



l -N(N-iy. 



(510) 



24 ti 2 

The deviation of the Hartree-Fock result from the exact result is a fraction 1/N of the 
energy and is small indeed in the large N limit, as expected from the validity condition 
(495)! 

There is a notable difference of translational properties however. Whereas the exact 
ground state (507) is invariant by a global translation of the positions of the particles, 
as it should be, the Hartree-Fock ansatz leads to condensate wavefunctions tp localized 
within the length / around some arbitrary point zq (around zq = in Eq.(491)): 

1 1 



il>zo(z) = 



(2/) 1 / 2 CQ8h[(z-Z )/l] 



with a spatial radius 



2h 2 



(511) 



(512) 



(N-l)mg ld 

The Hartree-Fock ansatz \N : ip) therefore breaks the translational symmetry of the 
system. 

5 The center of mass of the gas corresponds to a fictitious particle of wavevector K , where KK 
is the total momentum of the gas, and of position Z , where Z is the centroid of the gas. In the 
ground state |\&) the center of mass is completely delocalized with K = . The factor 1/L in |AT| 2 
originates from the normalization of the fictitious particle plane wave in the fictitious box of size L , 
(Z\K) = e lKZ fs/L . The more correct mathematical way (not used here) is to normalize in free space 
(no box) using the closure relation f dK\K)(K\ = Id , which amounts to replace L by 2tt . 
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Breaking a symmetry of the system costs energy, and this can be checked for the 
present translational symmetry breaking. As the center of mass coordinates Z, P of the 
N particles are decoupled from the relative coordinates of the particles we can write the 
total energy of the gas as the sum of the kinetic energy of the center of mass and an 
'internal' energy including the kinetic energy of the relative motion of the particles and 
the interaction energy. Whereas the exact ground state wavefunction has a vanishing 
center of mass kinetic energy, the symmetry breaking ansatz |iV : ^) contains a center 
of mass kinetic energy: 

E c ,. m . = (N :^\^\N (513) 

where mN is the total mass of the gas and P is the total momentum operator. Using 
the definition (497) , expanding the square of P , and using the fact that the soliton 
wavefunction ip has a vanishing mean momentum we obtain 

E c . .m. = (514) 

= y^ {n ~ 1)2 - (515) 

We see that E Cm0mmm accounts for half the energy difference between the exact ground state 
energy (506) and the Hartree-Fock energy (510). 

9.2.3 Physical advantage of the symmetry breaking description 

We now raise the question: is there a Bose-Einstein condensate in the one-dimensional 
free Bose gas with attractive interaction? To make things simple we assume that the 
gas is at zero temperature so that the N— particle wavefunction is known exactly, see 
Eq.(507). 

We start with a reasoning in terms of the one-body density operator (even if we know 
from the previous physical examples that this may be dangerous). Paraphrasing the 
usual three dimensional definition of a Bose-Einstein condensate in free space we put the 
one-dimensional gas in a fictitious box of size L and we calculate the mean number of 
particles n in the plane wave with vanishing momentum p = in the limit L — > +oo . 

The calculation with the exact ground state wavefunction has been done [85]. One 
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finds that n Q is going to zero as 1/L : 



2f> 2 

lim n Q L = C(N)— r . (516) 

l^+oo m\gi d \ 



The factor C(N) is given by 



N N (j - iv CAT - iV J r 1 i _1 

C(N) " ££ Nrlsdi S + 1 " t; " s'" + "1 (517) 

and converges to 7r 2 /2 in the large N limit, so that n no longer depends on TV in 
this limit. There is therefore no macroscopic population in the p = momentum state. 
One may then be tempted to conclude that there is no Bose-Einstein condensate, even at 
zero temperature, in the one-dimensional Bose gas with attractive contact interactions. 
However we have learned that a reasoning based on the one-body density matrix may 
miss crucial correlations between the particles, and that the symmetry breaking point of 
view may be illuminating in this respect. 

The translational symmetry breaking point of view approximates the state of the gas 
by the N -body density operator: 

^\N:^ )(N:^ \. (518) 

-L/2 L 

In the large N limit we expect this prescription to be valid for few-body observables. Of 
course for a N— body observable such as the kinetic energy of the center of mass of the 
gas, the results will be different, Eq.(515) for the symmetry breaking point of view vs. a 
vanishing value for the exact result. 

Let us test this expectation by calculating in the Hartree-Fock approximation the mean 
number of particles in the plane wave (z\k) — exp(ikz) / X 1 / 2 . Using the following action 
of the annihilation operator a k of a particle with wave vector k on the Hartree-Fock 
state: 

a k \N : ^o) = N l / 2 {k\^ Z0 )\N - 1 : (519) 

we obtain 

n* = N\{km 2 . (520) 

The momentum distribution of the particles in the gas in this approximation is simply 
proportional to the momentum distribution of a single particle in the solitonic wavefunc- 
tion ip ! It turns out that the Fourier transform of the 1 / cosh function can be calculated 
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exactly, and it is also a 1/ cosh function. We finally obtain: 



hf 1 k 2 % 2 1 



Uk "^m|^| C0S h 2 (f) (521) 

where / is the soliton size given in Eq.(512). For k = one recovers the large TV limit 
of the exact result (516). 

In more physical terms, one can imagine from Eq.(518) that a given experimental 
realization of the Bose gas corresponds to a condensate of N particles in the solitonic 
wavefunction (511), with a central position z being a random variable varying in an 
unpredictable way for any new realization of the experiment. There is therefore a Bose- 
Einstein condensate in the one-dimensional attractive Bose gas! 

An illustrative gedanken experiment would be to measure the positions along z of all 
the particles of the gas. In the symmetry breaking point of view the positions Zi, . . . ,z N 
obtained in a single measurement are randomly distributed according to the density 
\ipz \(z) = \tp(z — z Q )\ 2 where z Q varies from shot to shot as the relative phase of the two 
condensates did in the MIT interference experiment. As we know the exact ground state 
(507) we also know the exact N— body distribution function, \^(zi, . . . , z N )\ 2 . This is 
however not so easy to use! 

So we suggest instead to consider the mean spatial density of the particles knowing 
that the center of mass of the cloud has a position Z . In the exact formalism this gives 
[85]: 

P (z\z) = y - - - y rf^jvi^c^i, - - - , ^jv>r ^£ — ^,0^ ^ (-^ — ^ £ ( 522 ) 



2N I *z* (N-2)\ N\ , . fc/I . [ „ . 2N \z - Z\ 



I ^ (N -k-2)\(N + k)V ^v'-/— ^ v "'^iV-l I 

where I is the A -dependent length of the soliton (512), the integrals are taken in 
the range [-L/2, L/2] and L — »• +oo ; the factor L , compensating the one in the 
normalization factor of ^ , ensures that the integral of p(z\Z) over z is equal to A . 

In the symmetry breaking point of view the definition of p(z\Z) is similar to Eq.(522); 
the factor L cancels with the 1/L factor of Eq. (518). This leads to 

p sb (z\Z) = Jdz Jd Zl ...J dz N (f[ W** - *o)| 2 ) |£ S(z - s(z-^Y, 
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= N J dzi...J dz N \ip(zk)\ 2 ^ d^Z-z + zx-^J^Zn^ (523) 

where we have made the change of variables Zk — »• Zk + zq (which allows to integrate 
over z ) and we have replaced the sum over the indiscernible particles j by TV times 
the contribution of particle j — 1 . The multiple integral over the positions Zi, . . . ,z N 
can be turned into a single integral over a wavevector q by using the identity S(X) = 
J dq/(27r)exp(iqX) , allowing a numerical calculation of p sb (z\Z) . 

Does the approximate result (523) get close to the exact result for large N ? We 
compare numerically in figure 17 the exact density p(z\Z) to the symmetry breaking 
mean-field prediction p sb (z\Z) : modestly large values of N give already good agree- 
ment between the two densities. This validates the symmetry breaking approach for the 
considered gedanken experiment. 

What happens in the large N limit? In Eq.(523) each variable z k explores an interval 
of size rsj I so that the quantity (z\ + . . . + z N )/N has a standard deviation ~ l/y/N 
much smaller than / and can be neglected as compared to z\ inside the 6 distribution. 
This leads to 

p sb (z\Z)~N\iP Z0=z (z)\ 2 for y/N^>l (524) 

where the solitonic wavefunction ip zo =z is given in Eq.(511). Numerical calculation of 
p sb (z\Z) shows that Eq.(524) is a good approximation over the range \z — Z\ ~ I for 
N = 10 already! 
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Figure 17: For the ground state of the one-dimensional attractive Bose gas, position dependence 
of the mean density of particles knowing that the center-of-mass of the gas is in Z = . Solid 
line: exact result p(z\Z = 0) . Dashed line: mean- field approximation /? b (z\Z = 0) . The 
position z is expressed in units of the 'soliton' radius I given in Eq.(512), and the linear 
density in units of N/l . The number of particles is (a) N = 10 and (b) N = 45 . 
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